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I.     INTRODUCTION 

The  subject  of  this  dissertation  is  the  investigation  of  thermoacoustic  prime  movers  in 
an  annular  geometry.  What  sets  this  research  apart  from  previous  work  in  prime  movers  is 
the  nature  of  the  boundary  conditions.  To  our  knowledge  this  work  is  the  first  thorough 
investigation  of  thermoacoustic  prime  movers  with  periodic  boundary  conditions.  The 
primary  conclusion  drawn  from  this  work  is  that  a  full  understanding  of  the  eigenmodes  of 
the  system  are  required  to  design  a  functioning  annular  prime  mover. 

Thermoacoustic  heat  transport  is  a  process  through  which  an  acoustic  field  generates, 
or  is  generated  from,  a  flow  of  heat.  Thermoacoustic  engines  are  of  two  types:  heat  pumps 
and  prime  movers.  A  thermoacoustic  heat  pump  utilizes  a  standing  wave  to  transport  heat 
along  the  boundary  of  a  plate  situated  in  the  standing  wave.  The  acoustically  generated  heat 
flow  produces  a  temperature  gradient  across  the  stack.  In  other  words,  acoustic  energy  is 
converted  into  stored  thermal  energy.  In  contrast,  a  thermoacoustic  prime  mover  produces 
acoustic  work,  by  accepting  heat  from  a  high-temperature  source  and  transferring  it  to  a 
low-temperature  sink.  In  this  dissertation  a  prime  mover  in  an  annular  resonator  is 
investigated  both  theoretically  and  experimentally. 

To  illustrate  how  a  prime  mover  operates,  two  conventional  prime  mover 
configurations  are  described.  One  type  (as  is  shown  in  Fig.  1.1)  is  a  rigid-rigid  standing 
wave  tube  which  contains  a  stack  of  plates  called  the  prime  mover  stack  (or,  simply,  the 
stack),  which  is  in  thermal  contact  with  two  heat  exchangers.  In  operation,  a  temperature 
difference  is  imposed  between  the  two  heat  exchangers,  resulting  in  a  temperature  gradient 
along  the  stack.  When  the  temperature  gradient  in  the  stack  is  sufficiently  large,  the  gas  in 
the  resonator  oscillates  spontaneously  at  a  certain  frequency,  with  pressure  antinodes  at  the 
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closed  ends.  The  inception  of  spontaneous  oscillation  is  known  as  onset.  Once  onset  is 
reached,  the  acoustic  amplitude  in  the  resonator  grows  rapidly,  typically  reaching  several 
percent  of  the  mean  gas  pressure. 
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Figure  1.1  A  typical  rigid-rigid  prime  mover  configuration. 


The  second  prime  mover  configuration,  called  the  Hofler  tube,  is  a  rigid-open 
resonator  (as  is  shown  in  Fig.  1.2).  In  1983,  Hofler  designed  and  built  this  simple 
thermoacoustic  oscillator  to  demonstrate  some  thermoacoustic  phenomena  to  his  doctoral 
committee  at  the  University  of  California,  San  Diego.  The  sound  radiated  from  the  open 
end  of  the  pipe  is  impressively  loud,  about  100  dB  re  20  |iPa  a  meter  away  (Wheatley  et 
al.,  1985;  Swift,  1988). 
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Figure  1.2  The  Hofler  tube,  an  example  of  a  rigid-open  prime  mover. 


The  condition  necessary  for  a  prime  mover  to  reach  onset  is  that  the  amount  of  stored 
thermal  energy  converted  to  acoustic  energy  must  exceed  the  total  amount  of  acoustic 
energy  dissipated  by  thermal-viscous  losses  in  the  prime  mover.  This  ability  to  reach  onset 
is  determined  by  the  geometry  of  the  prime  mover,  in  particular,  the  position  of  the  stack 
relative  to  the  nearest  pressure  antinode.  The  hot  end  of  the  stack  is  closer  to  a  pressure 
antinode  than  is  the  ambient  end.  In  the  two  examples  of  conventional  prime  movers 
described  above,  the  node  and  antinode  positions  are  fixed  by  incorporating  well-defined 
acoustic  boundary  conditions  into  the  system.  This  is  an  important  common  feature  for  all 
conventional  prime  movers.  The  "built-in"  boundary  conditions  also  serve  as  convenient 
starting  points  for  computational  analysis  of  the  performance  of  a  prime  mover.  Because 
the  relative  positions  of  the  dominating  boundaries  and  the  stack  are  fixed,  the  stack 
position  is  optimized  for  only  one  acoustic  mode.  Also,  the  optimal  stack  position  for 
onset  is  not  necessarily  the  optimal  position  for  high  amplitude  performance.  One  of  the 
initial  motivations  for  studying  annular  prime  movers  was  to  see  if  it  would  self-optimize 
the  stack  location  relative  to  the  acoustic  field. 


Unlike  a  conventional  prime  mover,  an  annular  prime  mover  (as  is  shown  in  Fig.  1.3) 
does  not  have  the  typical  dominating  boundary  condition  to  force  a  pressure  or  velocity 
node  at  any  particular  position.  Design  and  analysis  of  a  functional  annular  prime  mover 
are  thus  made  more  difficult. 
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Figure  1.3  An  annular  prime  mover  configuration. 


In  a  uniform  cross-section  annular  resonator,  the  fundamental  longitudinal  mode 
corresponds  to  the  circumference  being  approximately  equal  to  one  acoustic  wavelength. 
There  are  two  degenerate  orthogonal  modes  that  satisfy  this  condition.  However,  these 
modes  have  no  preferred  orientation.  The  presence  of  a  nonuniformity  (for  example,  the 
stack)  breaks  the  degeneracy  resulting  in  a  frequency-splitting.  The  low  frequency  mode 
will  have  a  pressure  node  at  the  stack,  and  the  higher  frequency  mode  will  have  a  velocity 
node  at  the  stack.  (In  what  follows,  the  word  "frequency"  will  be  dropped  and  the  modes 


referred   to,    simply    as    "high"    and    "low".)       Neither    of    these    modes    supports 
thermoacoustics. 

Although  in  retrospect  the  answer  now  seems  obvious,  prior  to  a  thorough  study  of 
the  problem  one  might  have  argued  that,  once  it  becomes  very  strong,  the  thermoacoustic 
effect  may  dominate  the  situation  and  shift  the  orientation  of  the  acoustic  field  to  some 
preferable  position.  This  suggests  the  possibility  of  the  annular  prime  mover  achieving 
onset.  Even  if  the  uniform  cross  section,  single  stack  prime  mover  does  not  reach  onset, 
other  nonuniformities  in  cross  section  will  have  some  effect  on  the  spatial  distribution  of 
the  acoustic  field.  Hence,  by  placing  another  constriction  in  the  annulus,  the  previously 
mentioned  pressure  and  velocity  nodes  may  be  displaced  from  the  stack,  providing  the 
possibility  of  the  annular  prime  mover  reaching  onset.  These  are  some  of  the  concepts  that 
initiated  this  work. 

It  should  be  noted  that  a  rigid-rigid  prime  mover  can  be  considered  a  limiting  case  of  a 
constricted,  annular  prime  mover.  All  the  previous  work  on  thermoacoustic  prime  movers 
have  dealt  with  engines  that  have  some  sort  of  well-defined  boundary  conditions  imposed 
upon  them.  There  has  been  no  detailed  analysis  of  prime  movers  with  periodic  boundary 
conditions. 


A.  BACKGROUND 

Although  thermoacoustic  phenomena  have  been  observed  for  a  long  time,  significant 
advances  in  practical  thermoacoustics  are  relatively  recent.  The  earliest  example  of  a 
thermoacoustic  prime  mover  is  the  Sondhauss  tube  (Sondhauss,  1850).    Over  100  years 


ago,  it  was  observed  by  glassblowers  that  a  hot  glass  bulb  attached  to  a  cool  glass  tube 
sometimes  emitted  sound  at  the  tip  of  the  tube.  Sondhauss  quantitatively  investigated  the 
relationship  between  the  pitch  of  the  sound  emitted  and  the  dimensions  of  the  tube.  Lord 
Rayleigh  explained  the  Sondhauss  tube  quantitatively  in  1 896,  but  no  complete  theoretical 
analysis  of  these  phenomena  was  made  before  publication  of  the  series  of  papers  by 
Nikolaus  Rott  (Rott,  1969,  1980,  1983).  Later,  Wheatley  (1985)  and  others  at  Los 
Alamos  began  the  development  of  practical  thermoacoustics  devices. 


B.    PAST  RELEVANT  WORKS 

Past  works  relevant  to  this  research  fall  under  three  headings:   (1)  the  thermoacoustic 
prime  mover,  (2)  the  acoustic  Stirling  Engine,  and  (3)  the  annular  resonator. 

1.     The  Thermoacoustic  Prime  Mover 

The  potential  application  of  thermoacoustics  as  a  heat-driven  sound  source  has 
motivated  theoretical  developments  (Rott,  1969;  Wheatley  and  Cox,  1988;  Swift,  1988). 
Work  on  prime  movers  has  been  focused  on  either  rigid-rigid  or  rigid-open  prime  movers. 
Some  examples  are  outlined  below.  Miglori  and  Swift  (1988)  constructed  a  thermoacoustic 
prime  mover  that  used  liquid  sodium  as  its  working  fluid.  A  prime  mover  has  also  been 
used  as  a  heat-driven  sonar  projector  (Gabrielson,  1991).  Swift  (1992)  built  and  analyzed 
a  5-inch  thermoacoustic  prime  mover  which,  at  its  most  powerful  operating  point,  using 
13.8-bar  helium,  delivered  630  W  to  an  external  acoustic  load.  With  minor  modifications, 
this  device  was  used  by  Swift  (1994)  to  research  the  application  of  similitude  to  nonlinear 
thermoacoustics. 


The  Naval  Postgraduate  School  has  been  in  the  forefront  of  much  thermoacoustic 
research.  Atchley  and  the  author  first  built  and  tested  a  heat  driven  prime  mover  (Lin  and 
Atchley  1989),  which  generated  a  sound,  at  a  temperature  difference  of  453  °C,  with  a 
peak  acoustic  amplitude  of  7.9%  of  atmospheric  pressure.  Subsequently,  other  work  was 
done  on  this  device  (Atchley  1992,  1993;  Atchley  and  Kuo,  1994).  In  1993,  a 
thermoacoustic  prime  mover  was  constructed  by  Castro  and  Hofler  to  investigate  the 
performance  of  heat  exchangers  and  produced  peak  pressures  of  up  to  20%  of  the  mean 
pressure  (Castro  and  Hofler,  1993). 

DeltaE,  a  program  developed  at  Los  Alamos  National  Laboratory  by  Ward  and  Swift 
(1994,  1996)  has  been  applied  to  several  cases  in  the  design  and  analysis  of  thermoacoustic 
prime  movers.  DeltaE  stands  for  Design  Environment  for  Low  Amplitude  Thermoacoustic 
Engines.  It  solves  Rott's  wave  equation  for  a  geometry  defined  by  the  user.  For  example, 
Swift  (1992,  1996)  used  DeltaE  to  analyze  the  performance  of  a  5-inch  thermoacoustic 
engine.  Nessler  (1994)  used  DeltaE  to  compare  the  performance  of  the  pin  stack  with  a 
conventional  stack  in  a  thermoacoustic  prime  mover  (Swift  and  Keolian  1993).  Yang 
(1995)  and  Meng  (1996)  also  used  DeltaE  to  predict  onset  of  their  thermoacoustic  prime 
movers.  More  recently,  DeltaE  was  applied  to  analyze  the  efficiency  of  a  thermoacoustic 
prime  mover  with  a  pin  stack  (Gibson,  Nessler,  and  Keolian,  1997). 

2.     The  Acoustic  Stirling  Engine 

Ceperley  (1979,  1982  and  1985)  has  discussed  acoustic  engines  using  traveling 
waves.  As  recognized  by  Ceperley  and  several  other  authors  (Swift,  1988,  1995;  Hofler 
1988;  Organ  1992),  the  phasing  between  pressure  oscillations  and  velocity  within  a  Stirling 
engine  regenerator  is  the  same  as  those  of  an  acoustic  traveling  wave.  Ceperley  proposed 
the  replacement  of  the  function  of  the  pistons  in  a  Stirling  engine  by  acoustic  processes,  to 


form  what  he  called  a  traveling-wave  heat  engine.  One  of  his  concepts,  a  traveling-wave 
heat-driven  refrigerator,  is  shown  in  Fig.  1.4.  In  this  device,  one  regenerator  and  heat 
exchanger  set  functions  as  the  prime  mover,  adding  acoustic  power  onto  the  traveling  wave 
as  heat  flows  from  a  high-temperature  heat  source  to  a  room  temperature  heat  sink.  The 
other  set  functions  as  a  heat  pump,  using  acoustic  power  from  the  traveling  wave  to 
transport  heat  from  a  low  temperature  to  room  temperature.  His  work  was  partly 
responsible  for  motivating  the  author  to  embark  on  this  project.  However,  his  analysis  was 
overly  simplistic  and  he  never  built  a  working  device. 
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Figure  1.4     Ceperley's  design  for  a  traveling-wave  heat-driven  refrigerator.     The 
arrows  show  the  direction  of  wave  propagation. 

3.    The  Annular  Resonator 


Previous  work  on  annular  resonators  is  also  relevant  to  the  study  of  the  annular  prime 
mover.  The  eigenmodes  for  electromagnetic  wave  propagation  in  an  annular  cavity  are  of 


fundamental  interest  in  the  radio-frequency  (RF)  heating  of  tokamak  plasmas,  besides 
having  other  applications  in  the  microwave  circuit  theory  (Cap  and  Deutsch  1978,  1980; 
Janaki  and  Dasgupta,  1990;  Wu,  1992). 

Another  category  of  relevant  work  is  the  frequency  splitting  phenomenon  found  in 
different  fields.  As  early  as  1969,  Rudnick  et  al.  made  observations  of  a  superfluid  helium 
persistent  current  using  the  Doppler-shifted  splitting  of  an  azimuthal  resonant  fourth-sound 
mode  of  the  cylindrical  resonator  containing  liquid  helium  (Rudnick  et  al.  1969). 
Heiserman  investigated  the  persistent  currents  in  superleaks  in  contact  with  bulk  superfluid 
helium  using  the  Doppler  shifts  of  the  acoustic  modes  of  an  annular  resonator  partially 
packed  with  a  superleak,  and  with  a  simple  gyroscopic  technique  (Heiserman,  1975).  In 
plasma  physics,  in  the  presence  of  a  finite  plasma  current,  the  axisymmetric 
magnetoacoustic  wave  resonance  exhibits  a  frequency  splitting  for  a  finite  annular  mode 
number  between  oppositely  directed  traveling  waves  (Borg  and  Wit,  1991).  A  deformation 
of  the  cross  section  of  a  pinch  of  the  plasma  cross  section  also  induces  frequency  shift  and 
splitting  of  circular  modes  (Ring,  1988).  In  the  microwave  and  RF  field,  frequency 
splitting  can  also  be  found  in  a  nonuniform  microstrip  ring  resonator  (Wested  and 
Andersen,  1972),  in  an  asymmetrically  coupled  microstrip  ring  resonators  (Al-Charchafchi 
and  Dawson,  1990)  and  in  a  microstrip  rhombic  resonator  with  a  step  discontinuity  in  one 
of  the  arms  (Al-Charchafchi  and  Boulkos,  1990;  Al-Charchafchi  and  Schreck,  1994). 

C.    OUTLINE  OF  THIS  WORK 

This  work  constitutes  the  first  detailed  theoretical  and  experimental  investigation  of  a 
thermoacoustic  prime  mover  with  periodic  boundary  conditions.  There  are  five  significant 
aspects  to  this  research:     (1)  using  DeltaE  to  analyze  an  annular  prime  mover,  (2) 


developing  an  entirely  new  analysis  program  using  MATLAB,  (3)  designing,  building,  and 
experimentally  investigating  a  single  stack,  annular  prime  mover,  (4)  experimentally 
investigating  a  constricted,  single  stack  annular  prime  mover,  (5)  predicting  the 
performance  of  a  two  stack  annular  prime  mover. 

The  major  conclusions  are:  (1)  A  single  stack  annular  prime  mover  will  not  reach 
onset  because  the  eigenmodes  of  the  system  do  not  support  thermoacoustic  growth.  (2)  A 
constricted  annular  prime  mover  will  reach  onset  because  the  constriction  forces  dominating 
boundary  conditions  that  alter  the  eigenmodes.  (3)  A  two  stack  prime  mover  is  predicted  to 
reach  onset  because  one  of  the  eigenmodes  of  the  symmetric  system  does  support 
thermoacoustic  growth. 
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II.    THEORY 

This  chapter  is  divided  into  four  parts.  The  first  part  is  a  brief  review  of  the  aspects  of 
thermoacoustics  pertinent  to  the  annular  prime  mover  problem.  The  review  draws  heavily 
from  Swift's  work  (Swift,  1988  and  1997),  and  is  intended  to  introduce  the  concepts 
required  to  understand  the  techniques  used  to  analyze  the  annular  prime  mover.  One  of  the 
central  questions  to  be  addressed  in  this  dissertation  is  the  ability  of  an  annular  prime  mover 
to  reach  onset.  The  dependence  of  quality  factor  Q  of  the  prime  mover  on  the  temperature 

difference  AT  applied  to  the  stack  is  a  key  measure  of  this  ability.    In  this  research,  the 

predicted  Q  will  be  determined  from  the  complex  eigenfrequency  of  the  prime  mover. 
Therefore,  following  the  review  of  thermoacoustics,  the  complex  eigenfrequency  will  be 
introduced. 

Another  important  aspect  of  this  dissertation  is  the  numerical  analysis  of  the  annular 
prime  mover.  The  techniques  used  will  be  discussed  in  the  last  part  of  this  chapter.  Some 
basic  properties  of  the  two-point  boundary-value  problem  are  discussed,  followed  by  a 
discussion  of  the  advantages  and  disadvantages  of  using  DeltaE  with  the  annular  prime 
mover  geometry.  A  MATLAB  numerical  analysis  program  is  then  described  that  is  more 
tailored  to  the  annular  prime  mover  than  is  DeltaE.  The  validation  of  this  program  is 
discussed  and  the  application  of  the  MATLAB  program  to  our  annular  prime  mover  is 
presented. 

This  chapter  concludes  with  a  discussion  of  the  end  corrections  for  constrictions  in  the 
cross  section  of  annular  resonators. 
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A.    REVIEW  OF  THERMOACOUSTICS 

Thermoacoustics  is  essentially  the  study  of  the  acoustics  of  a  fluid  in  which  there 
exists  a  temperature  gradient.  The  acoustics  of  such  a  fluid  is  contained  in  three  differential 
equations  which  relate  the  complex  acoustic  pressure  and  volume  velocity,  and  the 
temperature  and  enthalpy  of  the  fluid.  This  method  was  first  described  by  Wheatley  et  al. 
(Wheatley  et  al.,  1983).  The  derivations  of  these  relationships  are  summarized  in  this 
section.  Once  these  concepts  have  been  introduced,  a  description  of  the  basic  scheme  used 
to  analyze  the  performance  of  an  annular  prime  mover  is  given. 

This  section  is  essentially  a  summary  of  pertinent  aspects  of  thermoacoustics  given  in 
the  articles  by  G.  W.  Swift  (Swift,  1988  and  1997).  Swift's  quantitative  thermoacoustic 
analysis  is  based  on  Rott's  wave  equation  and  energy-flux  equation,  which  result  from  the 
momentum,  mass  continuity,  and  energy  equations  in  the  acoustic  approximation.  The 
results  are  valid  for  arbitrary  phasing  between  the  acoustic  pressure  and  velocity;  in  other 
words,  for  both  standing-wave  and  traveling-wave  systems.  The  results  presented  in  this 
section  serve  as  the  basis  for  the  numerical  analysis  used  in  this  dissertation. 

First,  the  notation  that  will  be  used  throughout  this  discussion  is  established.  The 
acoustic  wave  is  considered  to  propagate  in  the  x  direction  within  a  duct  of  constant  cross 
section.  The  usual  complex  notation  is  adopted  for  time-oscillatory  acoustic  quantities. 
Also,  expansions  to  first  order  in  the  acoustic  variables  is  assumed  to  suffice.  Thus,  the 
pressure,  velocity  and  temperature  can  be  written 

p(x)=pm+P](x)c"*,  (2.1) 

u(x,  y,  z)  =  ux  (x,  y,  z)  e'"  ,  (2.2) 

and 
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T(x,  ,y,  z)  =  T  m  (x)  +  T,  (x,  y,  z)  eia" ,  (2.3) 

where  i  =  V-l  and  w,  is  the  x-component  of  v .  The  mean  values  (subscript  m)  are  real, 
but  the  acoustic  variables  (subscript  1)  are,  in  general  complex,  to  account  for  the  relative 
phasing  of  those  oscillating  quantities.  It  is  assumed  that  the  mean  fluid  velocity  um  =  0. 
All  the  time  dependence  appears  in  the  emt  term,  with  co  =  2  n  f  being  the  angular 
frequency.  For  clarity,  we  consider  here  only  the  case  of  large  solid  heat  capacity,  so  that 
the  temperature  of  the  solid  material  in  the  stack  is  simply  Tm(x),  independent  of  t,  y,  and  z. 
The  conservation  of  momentum  of  an  incompressible  viscous  fluid  is  ( Landau,  1975) 


^  +  (v  .grad)v  =    -  1  grad/?  +^  V2v,  (2.4) 


where  v  is  the  velocity,  p  is  the  density,  p  is  the  pressure  and  fi  is  the  dynamic  viscosity  of 
the  fluid  in  the  resonator.  This  is  called  the  Navier-Stokes  equation.  To  develop  a 
quantitative  understanding  including  viscous  effects,  we  begin  by  finding  the  y  and  z 
dependence  of  the  velocity.  Because  of  the  viscous  boundary  layer  in  the  y  and  z  direction 
all  other  viscous  derivatives  can  be  neglected  compared  to  jid2u^dy2  and  fid2uxldz2. 
Equation  (2.4)  then  reduces  to 


i(opmux  =    -— -  +  p. 
ax 


dp]  fd2Uj         d2Uj 


ydy2  dz2  j 


(2.5) 


Equation  (2.5)  is  an  ordinary  differential  equation  for  ux{y,  z).  The  boundary 
condition  at  the  solid  surface  is  «,  =  0.  With  this  boundary  condition  Eq.  (2.5)  can  be 
solved  to  yield 
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ux(x,  y,  z)=   —  [  l-hv(y,  z)]-^-  ,  (2.6) 

copm  dx 

where  hv(y,  z)  depends  on  the  specific  geometry  under  consideration.  The  volume  velocity 
is  the  spatial  average  over  v  and  z  of  Eq.  (2.6)  and  is  thus 


Urix)  =  i^L(l-fv)*EL,  (2.7) 

copm  dx 


where  /v  is  the  spatial  average  of  hv(y,  z).    For  the  case  of  parallel-plate  geometry  of  the 
stack  used  in  this  work,  it  can  be  shown  that  (Swift  1988;  Arnott  et  al.  1991) 


,    ,  x        cosh[(l  +  /)y/<5v] 

K(y)=         "     - ,:.,  (2.8) 

cosh[(l  +  0y0/<5j 


and 


tanh[(7  +  i)y0/gy] 


where  <SV  =  A/2iu/pmO)  is  the  fluid's  viscous  penetration  depth  and  v0  is  plate  half-gap. 

Calculation  of  the  oscillating  fluid  temperature  is  accomplished  using  only  the  first 
order  terms  in  the  general  equation  of  heat  transfer  (Landau,  1975).  The  entropy  is 
expressed  in  terms  of  the  acoustic  pressure  and  temperature  of  the  fluid  and  the  resulting 
ordinary  differential  equation  is  solved  subject  to  the  appropriate  boundary  conditions.  In 
the  presence  of  a  temperature  gradient  along  the  duct,  the  solution  is  (Swift  1988  and  1997) 


T](x,y)=    -L-[l-hK]Pl  +       — j—      ^Ml-V  °K   )Ux,      (2.10) 
Pmcp  <w(l-/v)Agas   dx  1  -  a 


where  (J  =  cp  jh/k  is  the  Prandtl  number,  and  cp  is  the  isobaric  heat  capacity  per  unit  mass. 
The  spatial  average  in  (v,  zj-direction  of  Eq.  (2.10)  can  be  written 
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Ti(x)=  -l~[l-fK]P]  +  l  ^k(1.A  -o-/v 


where /K  is  the  spatial  average  of  hK.  The  functions  hK  and/K  are  the  same  as  /zv  and/v,  only 


with  <5V  replaced  by  SK=J2K/pmcpa>,  the  thermal  penetration  depth,  where  /cis  the  thermal 
conductivity  of  the  fluid.  The  function  fx  and/v  are  very  important  in  thermoacoustics  and 
describe  how  viscous  and  thermal  processes  affect  the  oscillatory  velocity  and  temperature 
in  the  acoustic  field.  A  graph  of/v  as  a  function  of  y/5v  is  shown  in  Fig.  2.1. 
The  equation  of  continuity  is  (Landau,  1975) 


-£  +  V  •  (pv)  =  0.  (2.12) 

at 

To  first  order,  the  x-component  of  Eq.  (2. 12)  can  be  reduced  to  the  form 


.       d(piio=o  (2i3) 

dx 
Equation  (2. 13)  can  be  rewritten  in  terms  of  volume  velocity 


1    d(pmJ7,) 

ifflp,  +  -r —j "   =  0-  (2-14) 

Combining  Eqs.  (2.10)  and  (2.13)  yields 


fi   =    -   ^[l  +  {Y-mPl  +  J^Lj-^,  ,  (2,5) 

dx  pmc2  (l-fv){l-o)Tm    dx 


where  c  is  the  speed  of  sound  in  the  gas,  and  7  is  the  ratio  of  isobaric  to  isochoric  specific 
heats. 
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We  pause  to  discuss  the  results  derived  in  the  last  several  paragraphs.  Equation  (2.7) 
shows  that  the  gas  which  is  much  farther  than  <5V  from  the  nearest  solid  surface  experiences 
essentially  no  viscous  shear  and  moves  with  a  velocity  that  depends  only  on  x.  The  gas 
that  is  much  closer  than  <5V  is  nearly  at  rest.  Gas  that  is  between  these  extremes  moves  with 
a  reduced  velocity  amplitude  and  significant  phase  change.  The  terms  in  Eq.  (2.11) 
indicate  that  there  are  two  contributions  to  the  oscillatory  temperature.  The  first  term  is 
simply  due  to  the  adiabatic  compressions  and  expansions  of  the  fluid.  The  second  term 
comes  from  the  convection  of  gas  parcels  along  a  temperature  gradient.  The  net 
temperature  oscillation  is  just  a  linear  combination  of,  these  two  effects.  It  should  be 
pointed  out  that/K  has  the  same  functional  dependence  as  y/8K  as  is  portrayed  in  Fig.  2.1 
for/. 

Equation  (2.15)  has  an  easy  physical  interpretation.  It  indicates  that  dU\/dx  comes 
from  a  density  oscillation  in  the  fluid.  This  density  change  can  be  caused  by  a  pressure 
change  or  by  convection  of  gas  along  the  temperature  gradient. 
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Figure  2.1   The  real  and  imaginary  parts  and  magnitude  of/v  as  functions  of  the  ratio 
y/£v  for  the  case  of  parallel  plate  geometry. 

Next,  following  the  derivation  by  Swift  (1988  and  1997),  the  time  averaged  energy 
flux  along  the  stack  can  be  found  as 


1  „ 
H     =   -  Re 

2         2 


-r*(x     a-/;) 


\ 


(1  +  (T)(1-/V) 


v  >J 


2fi)A^(l-<72)|l-/v| 


Im(/.  +  O//)^-  [\as  K  +  Km  Ksoltd)d 


(2.16) 


dx 


dx 


where  the  superscripts  *  represent  the  complex  conjugate  of  a  complex  quantity,  and  KsoVld 
is  the  thermal  conductivity  of  the  solid.  On  the  basis  of  Swift's  assumptions,  in  steady 
state  for  an  engine  without  lateral  heat  flows  to  the  surroundings,  H2  along  x  must  be  a 
constant  throughout  the  stack. 
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To  numerically  analyze  the  performance  of  thermoacoustic  engines,  it  is  recognized 
that  Eqs.  (2.7),  (2.15)  and  (2.16)  comprise  a  set  of  three  second-order  coupled  complex 
differential  equations  in  the  five  variables  Re(/?,),  Im^,),  Re(£/i),  Im(£/i)  and  T  .  These 
equations  are  the  basis  of  both  DeltaE  and  the  MATLAB  program  presented  later  in  this 
chapter. 

In  order  to  obtain  a  more  useful  form  of  the  differential  equation  for  p,(x),  Eq.  (2.7) 
can  be  written  as 


*■         i(°P"        Utix).  (2.17) 


dx        A(\-fv) 


In  a  duct  region,  the  integration  is  governed  by  Eqs.  (2.15)  and  (2.17)  and  the  temperature 
gradient  dTJdx  is  determined  by  the  measured  temperature  profile.  It  is  assumed  that  there 
is  no  temperature  gradient  across  the  heat  exchanger.  As  a  result,  in  a  heat  exchanger, 
integration  is  carried  out  by  Eq.  (2.17)  and  a  simplification  of  Eq.  (2.15) 

f-   =    -  ^%[l  +  (r-l)A]p,  .  (2.18) 

dx  pmc 

As  for  the  stack  region,  integration  is  based  on  Eqs.  (2.15),  (2.17),  and  (2.19).  In 
other  words,  the  time  averaged  energy  flux  H2  is  only  used  in  the  stack  region  and  is  a 
constant  throughout  the  stack  region. 

The  MATLAB  numerical  analysis  method  used  in  this  work  starts  with  specifying 
Re(/?,),  Im(p,),  Re((/i),  lm(U\),  Tm,  Re(/),  and  Im(/)  at  a  point  in  the  duct.  The  complex 
pressure  and  volume  velocity  just  outside  the  starting  end  of  the  ambient  heat  exchanger  can 
be  determined  by  integrating  Eqs.  (2.15)  and  (2.17)  in  the  boundary  layer  approximation 
limit.   Next,  the  complex  pressure  and  volume  velocity  can  be  determined  just  inside  the 
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starting  end  of  the  ambient  heat  exchanger  by  invoking  continuity  of  pressure  and  volume 
velocity  at  the  end  of  the  ambient  heat  exchanger.  Then,  Eqs.  (2.17)  and  (2.18)  are  used  to 
integrate  over  the  ambient  heat  exchanger.  Continuity  of  pressure  and  volume  velocity  is 
used  to  enter  the  stack.  A  value  of  H2  is  specified  to  get  drm/dx  at  the  beginning  of  the 
stack.  The  solution  in  integrated  to  the  heat  exchanger.  Again  Eqs.  (2.17)  and  (2.18)  are 
integrated  over  the  hot  heat  exchanger.  Finally,  the  pressure,  volume  velocity,  and  the 
temperature  at  the  exit  of  the  hot  heat  exchanger  are  matched  to  those  at  the  starting  point. 

B.     COMPLEX  EIGENFREQUENCY 

The  transition  to  onset  in  a  prime  mover  is  conveniently  discussed  in  terms  of  the 
quality  factor  Q  which  can  be  defined  as  the  ratio  of  2n  times  the  energy  stored  in  the 
resonator  to  the  energy  dissipated  per  acoustic  cycle.  Letting  W  denote  the  dissipated 
power, 

Q  =  ^r2--  (2.19) 

W 

Both  W  and  Esl  are  second  order  in  the  acoustic  amplitude.  A  typical  prime  mover  is 
comprised  of  five  sections:  the  ambient  duct,  the  ambient  heat  exchanger,  the  prime  mover 
stack,  the  hot  heat  exchanger,  and  the  hot  duct.  W  is  then  the  sum  of  the  power  dissipated 
in  the  five  individual  sections  (Lin,  1989;  Atchley,  1992  and  1994) 

W   =    Wamb+   Wambhx    +    Wpms    +    Wnothx    +    Who,.  (2.20) 

The  subscripts  amb,  amb  hx,  pms,  hot  hx,  and  hot  refer  to  the  ambient  duct,  the  ambient 
heat  exchanger,  the  prime  mover  stack,  the  hot  heat  exchanger,  and  the  hot   duct, 
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respectively.  The  reader  should  refer  to  Atchley  (Atchley  et  al.,  1992)  and  Lin  (Lin,  1989) 
for  a  full  explanation. 

Instead  of  Q,  we  prefer  to  use  the  reciprocal  of  Q,  because  l/Q  converges  to  zero  at 
onset  rather  than  diverging  to  infinity.  Also,  it  is  this  quantity  that  is  proportional  to  the 
acoustic  power  output  of  the  prime  mover.  It  is  worth  mentioning  that  l/Q  is  a  function  of 
the  prime  mover  geometry,  the  thermophysical  properties  of  the  gas,  pm,  and  the 
temperature  difference.  As  the  temperature  difference  along  the  stack  increases  from  zero, 
the  thermal  losses  in  the  stack  decrease  and  eventually  become  negative,  representing  gain. 
In  other  words,  below  onset  l/Q  is  positive,  at  onset  l/Q  is  zero  and  above  onset  l/Q  is 
negative.  A  complete  evolution  of  Q  through  onset  for  a  prime  mover  is  discussed  by 
Atchley  (Atchley,  1994). 

In  this  work,  the  Q  is  determined  by  the  complex  eigenfrequency.  Complex 
eigenvalues  occur  when  a  system  has  underdamped  modes.  The  excitation  of  a  resonator 
evolves  in  time  as  exp(/<w0f)  where  the  eigenfrequency  is  (Swift,  1988;  Gamaletsos  1993; 
Arnott  et  al.,  1994) 


co0  =  2Kf0  +  i^-,  (2.21) 

where  f0  is  the  real  resonance  frequency.  The  quality  factor  Q  thus  can  be  found  by 


Re(a)0 ) 

Q  =  ttV\  •  (2-22) 

2  Im(ft)0 ) 


It  should  be  noted  that  the  complex  eigenfrequency   approach   accounts   for  power 
generation  in  the  stack  and  dissipation  everywhere  in  the  resonator.    If  this  approach  were 


20 


not  used,  an  artificial  acoustic  source  would  have  to  be  introduced  into  the  system  to  deliver 
power  to  match  boundary  conditions  (Choe,  1997). 

C.    THE  ANNULAR  GEOMETRY  FOR  A  PRIME  MOVER 

As  discussed  in  the  previous  chapter,  the  major  difference  between  conventional  prime 
mover  geometry  and  that  of  an  annular  prime  mover  is  that  there  are  no  typical  dominating 
boundary  conditions  in  the  latter.  Instead,  an  annular  prime  mover  satisfies  the  periodic 
conditions 

px{r,  6)  =  P](r,  0  +  2  n),  (2.23) 

and 

U\  (r,  6)  =  Uiir,  0  +  2  it).  (2.24) 

At  first  glance,  this  difference  does  not  seem  to  cause  much  difficulty  in  the  design  and 
analysis  of  an  annular  prime  mover,  until  one  realizes  the  fact  that  almost  all  the  previous 
work  on  conventional  prime  movers  employ  some  dominating  boundary  condition  as  a 
convenient  starting  point  of  their  analysis.  Elimination  of  these  boundary  condition  makes 
research  on  an  annular  prime  mover  considerably  different  from  that  of  a  conventional 
prime  mover. 

It  can  be  shown  that  an  annular  resonator  can  be  modeled  as  a  straight  duct  having  an 
effect  length  Leff  that  is  related  to  the  inner  and  outer  radii  of  the  annulus  and  the  particular 
mode  under  consideration  (Choe,  1997).  Therefore,  we  will  treat  the  annular  prime  mover 
as  a  straight  duct  having  periodic  boundary  conditions.  L^can  be  obtained  by  first  writing 
the  general  solution  of  the  wave  equation  in  cylindrical  coordinates  which  includes  a  linear 
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combination  of  the  cylindrical  Bessel  and  Neumann  functions.  Note  that  the  Neumann 
function  must  be  retained  because  the  origin  is  excluded  in  the  annular  resonator.  The 
effective  circumference  of  the  resonator  can  be  then  obtained  by  applying  rigid  boundary 
conditions  at  the  inner  and  outer  walls  of  the  annulus.  The  annulus  used  in  the  experiments 
presented  later  has  inner  and  outer  radii  of  10.0  and  15.3  cm,  respectively.  The  height  of 
the  duct  is  5.0  cm.  The  cutoff  frequency  for  cross  modes  is  approximately  6.6  kHz.  We 
are  interested  in  the  fundamental  longitudinal  mode  (A  =  Leff)  corresponding  to  frequencies 
in  the  400  ~  450  Hz  range,  well  below  cutoff. 

D.    COMPUTER  SIMULATION  OF  THE  ANNULAR  PRIME  MOVER 

This  subsection  discusses  the  numerical  analysis  techniques  for  the  annular  prime 
mover.  For  the  purpose  of  numerical  integration  in  a  thermoacoustic  engine,  Eqs.  (2.7), 
(2.15),  and  (2.16)  can  be  rewritten  as  a  set  of  five,  first  order,  ordinary  differential 
equations  in  terms  of  five  independent  acoustic  variables:  Re(/?,),  Im(/?,),  Re(f/i),  Im(£/i) 
and  7L. 


dRe(Ui) 


dx 


d*e{Pi) 
dx 

d  Im(p, ) 


dx 


=  Re 


=  Im 


icop  A       — 
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(2.26) 


[l  +  (y-l)/Jp,  +    ,      \ ;       : -Ui\,      (2.27) 


22 


dhn(Ui) 


dx 


and 


dT 

m 

dx 


-  ^%+(r-i)/,]A  + 


U-/v) 


lffi& 


(i-/v)(i-o-)rM  dx 


//; 


Agas      2Agas 


Re 


P.  tf 


1- 


(fK  -ft)     \ 


P    c    \U\\ 

r  m      p   I      '  I 


(l  +  cr)(l-/;) 

Ma  +  <r/v*) 


2Ae2a5(l-cJ2)|l-/v 


g<" 


Re(ft)) 


4 


io/i'd 


(2.28) 


(2.29) 


Our  intention  is  to  solve  this  set  of  coupled  ordinary  differential  equations  for  Re(pj), 

Im(p,),  Re(  U\),  Im(  f/i)  and  Tm,  subject  to  the  periodic  boundary  conditions  at  the  ends  of 
the  prime  mover.  This  is  the  same  basic  method  employed  by  DeltaE. 

1.    Problem  Statement 

When  ordinary  differential  equations  are  required  to  satisfy  boundary  conditions  at 
more  than  one  value  of  the  independent  variable,  the  resulting  problem  is  called  a  two  point 
boundary  value  problem.  Determining  the  deflection  of  a  bar  rigidly  fixed  at  both  ends  is  a 
typical  example  where  the  conditions  specified  are  the  deflections  of  the  elastic  curve  at  the 
supports.  Heat-flow  problems  often  fall  into  this  class  when  temperature  or  temperature 
gradients  are  given  at  two  points.  The  problem  encountered  here  is  of  this  kind,  too. 

A  major  distinction  between  initial  value  problems  and  two  point  boundary  value 
problems,  as  is  pointed  out  by  Press  (Press  et  al.,  1992),  is  that  in  the  former  case  one  is 
able  to  start  an  acceptable  solution  at  its  beginning  (initial  values)  and  just  march  it  along  by 
numerical  integration  to  the  end  (final  values).  However  in  the  latter  case,  the  boundary 
conditions  at  the  beginning  point  do  not  determine  an  unique  solution  to  start  with.  For  this 
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reason,  two  point  boundary  value  problems  require  considerably  more  effort  to  solve  than 
do  initial  value  problems.  Keller  (Keller,  1992)  has  discussed  the  existence  and 
uniqueness  theory  for  two  point  boundary  value  problems. 

Two  different  numerical  approaches  to  the  annular  prime  mover  are  applied  in  this 
dissertation:  DeltaE  and  a  MATLAB  program. 

2.    DeltaE 

DeltaE  is  a  computer  program,  developed  by  Bill  Ward  and  Greg  Swift  at  Los  Almos 
National  Laboratory,  for  modeling  and  designing  thermoacoustic  and  other  one- 
dimensional  acoustic  apparatus.  Basically,  DeltaE  solves  the  one-dimensional  wave 
equation  based  on  the  usual  low-amplitude  acoustic  approximation.  It  numerically 
integrates  the  wave  equation  in  a  gas  or  fluid,  in  a  geometry  provided  by  the  user  as  a 
sequence  of  segments,  such  as  ducts,  transducers,  compliances,  heat  exchangers  and 
stacks.  It  uses  continuity  of  oscillating  pressure,  oscillating  volume  velocity,  and  mean 
temperature  to  pass  from  the  end  of  one  segment  to  the  beginning  of  the  next.  It  uses  the 
appropriate  wave  equation  and  temperature  equation  for  each  segment.  The  iteration  is 
controlled  by  global  parameters  such  as  frequency  and  mean  pressure,  and  by  local 
parameters  such  as  the  geometry  of  a  segment  and  enthalpy  flow.  By  defining  a  series  of 
inputs,  global  parameters  (the  guess  vector),  the  desired  outputs,  and  boundary  conditions 
(the  target  vector),  the  problem  is  solved  iteratively.  Guesses  are  updated  until  targets  are 
achieved.  The  number  of  elements  in  the  target  vector  must  equal  the  number  of  elements 
in  the  guess  vector,  otherwise  the  system  is  over-  or  under-determined.  When  the 
computed  and  specified  targets  agree  to  within  a  specific  tolerance,  the  solution  is  said  to 
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converge.  All  calculations  in  DeltaE  are  performed  in  double  precision  and  the  user  may  set 
the  error  tolerance  for  matching  targets. 

In  general,  because  a  pass  of  DeltaE' s  integration  does  not  reach  the  end  with  targeted 
values  of  all  variables,  a  shooting  method  is  used  to  adjust  chosen  initial  variables  in  order 
to  hit  the  target  values.  The  user  provides  guesses  for  chosen  initial  variables.  A 
successful  convergence  of  DeltaE  depends  heavily  upon  a  good  choice  of  guess  and  target 
members.  A  good  guide  to  the  choice  of  guess  and  target  variables  is  discussed  by  Ward 
and  Swift  (Ward  and  Swift,  1996). 

DeltaE  is  very  versatile  due  to  the  fact  that  the  elements  of  the  guess  vector  are  not 
limited  to  the  conventional  choices  consisting  of  real  and  imaginary  parts  of  pl  and  U\. 
Any  variables  that  have  an  effect  on  the  target  vector  variables  can  be  used  as  a  guess.  This 
feature  enables  DeltaE  to  compute  resonance  frequency,  temperature,  and  geometrical 
dimensions  in  order  to  satisfy  specified  boundary  conditions.  DeltaE  is  very  efficient  at 
converging  to  solutions  for  some  complicated  acoustic  systems,  however,  it  knows  nothing 
about  acoustics  or  physics.  For  example,  it  does  not  recognize  that  negative  frequencies, 
negative  pressures  or  negative  lengths  are  physically  meaningless.  It  simply  does  the 
integration.  Therefore,  the  reasonableness  of  the  solutions  produced  will  almost  always 
depend  on  the  quality  of  the  initial  guess  vector  and  choices  of  the  guess-target  vector. 

To  apply  DeltaE  to  the  annular  prime  mover  the  boundary  conditions  are  incorporated 
into  DeltaE' s  target  vector.  The  unknown  conditions  at  the  BEGIN  segment,  which  DeltaE 
is  supposed  to  find,  are  in  DeltaE' s  guess  vector.  To  implement  the  periodic  boundary 
conditions  into  DeltaE,  a  SOFTEND  segment  is  inserted  right  after  the  initial  BEGIN 
segment  in  the  input  file  of  the  DeltaE  program.  Another  SOFTEND  segment  (or  they  can 
both  be  HARDEND  segments)  is  used  at  the  end  of  the  model.  Four  DIFFTARGETs  are 
specified  which  match  the  amplitude  and  phase  of  px  and  U\  at  these  two  SOFTEND 
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segments  (Ward,  1997).  The  SOFTEND  segments  serve  only  as  a  tool  to  calculate  the 
outputs  and  will  be  ignored  in  the  target  vector.  SOFTENDs  do  not  force  the  impedance  to 
be  zero.  A  given  temperature  across  the  prime  mover  stack  is  achieved  by  selecting  the 
temperature  at  the  hot  heat  exchanger  as  a  target  and  the  amount  of  heat  input  to  the  ambient 
heat  exchanger  as  a  guess  vector.  We  now  have  five  targets  and  only  one  guess;  four  more 
guesses  are  still  required.  The  other  four  guesses  have  been  selected  to  be  the  amplitude 
and  phase  of  p,  and  U\  at  the  beginning  segment.  These  guesses  and  targets  (as  are  shown 
in  Appendix  B)  are  summarized  in  Table  2.1. 


Guess   vector 

Target  vector 

1.  Ip,  (begin)l 

1.  Ip,(begin)l  -  lp,(end)l  <  e 

2.  Phase[p,  (begin)] 

2.  Phase[p[(begin)]  -  phase[p,(end)]  <  £ 

3.  \U\  (begin)l 

3.  lJ7i  (begin)!  - 1 Z7 1  (end)l  <  e 

4.  Phase[£/i  (begin)] 

4.  Phase[  U\  (begin)]  -  phase[  U\  (end)]  <  e 

5.  Heat  input  at  ambient  heat  exchanger 

5.  Hot  heat  exchanger  temperature 

Table  2.1  Summary  of  the  guess  and  target  vectors  in  the  DeltaE  input  file  for  the  annular 
prime  mover,  e  represents  the  tolerance  for  convergence. 

The  important  parameters  for  this  work  are  the  Q,  the  resonance  frequency  and  the 
mode  shape  of  the  prime  mover.  To  attain  these  desired  quantities,  more  deliberation  has  to 
be  made  in  determining  the  input  file.  Two  methods  can  be  utilized  to  find  the  Q  with 
DeltaE  (Ward,  1997).  In  the  first  method,  the  Q  can  be  obtained  indirectly  by  adding  an 
artificial  driver  into  the  system,  sweeping  its  frequency,  and  computing  the  frequency 
response  of  the  system.    The  resonance  frequency  and  half  power  points  can  be  extracted 
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from  this  analysis.  An  alternate  method  is  to  compute  the  ratio  of  the  stored  energy  to  the 
dissipated  energy.  The  stored  energy  is  computed  by  adding  an  artificial  driver, 
determining  the  pressure  and  velocity  throughout  the  system  and  integrating  the  energy 
density.  The  dissipated  energy  is  computed  from  the  work  done  by  the  driver  to  excite  the 
system.  Since  DeltaE  only  computes  the  pressure  and  velocity  at  the  exit  of  a  segment,  the 
system  must  be  broken  into  many  segments  to  get  a  reasonable  approximation  to  the  stored 
energy.  Thus  the  frequency  response  method  was  chosen  to  find  the  Q  and  resonance 
frequency.  A  side  branch  electroacoustic  transducer  was  used  as  a  driver,  which  has 
frequency  independent  parameters.  A  MATHEMATICA  program  (refer  to  Appendix  A) 
was  also  developed  to  compute  a  least  squares  fit  to  the  frequency  response  to  find  Q  and 
resonance  frequency.  After  finding  resonance  frequency,  the  system  is  then  driven  at 
resonance.  To  attain  a  detailed  modal  shape,  the  prime  mover  has  to  be  divided  into  many 
small  segments.  The  modal  shape  is  a  plot  of  the  pressure  at  the  end  of  each  segments  of 
the  converged  solution  as  a  function  of  segment  location. 

DeltaE  is  a  very  useful  tool  to  predict  how  a  given  thermoacoustic  apparatus  will 
perform,  or  for  helping  the  user  design  an  apparatus  to  achieve  a  desired  performance. 
However,  it  does  have  several  disadvantages  when  applied  to  the  annular  problem.  As  was 
mentioned  in  Chapter  I,  in  the  presence  of  a  stack,  an  annular  prime  mover  exhibits 
frequency  splitting.  The  degree  to  which  a  particular  mode  is  excited  depends  on  the 
location  of  the  driver  relative  to  the  stack.  Because  of  dissipation,  Q  for  each  mode  is 
finite,  which  means  that  the  two  modes  may  overlap  in  the  frequency  domain.  In  terms  of 
a  driven  acoustic  system,  this  manifests  itself  as  a  multimode  excitation.  When  only  one 
mode  is  excited,  it  is  straightforward  to  find  Q  and  resonance  frequency  by  a  least  squares 
fit  to  the  frequency  response  curve.    However,  in  the  case  when  the  driver  excites  two 
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mode  simultaneously,  some  other  technique,  for  example  pole-zero  analysis,  must  be  used 
to  determine  Q  and  resonance  frequency  for  each  individual  mode. 

Another  disadvantage  comes  from  trying  to  map  out  the  mode  shape  with  DeltaE.  In 
DeltaE,  values  of  acoustic  pressure  are  only  provided  at  the  end  of  each  segment. 
Therefore,  to  obtain  a  detailed  mode  shape,  the  prime  mover  must  be  divided  into  many 
small  segments.  This  inherent  property  makes  it  is  laborious  and  impractical  to  use  DeltaE 
to  find  the  mode  shape. 

The  other  disadvantage  is  that  currently  there  is  no  straightforward  way  to  apply  an 
externally-imposed  temperature  gradient  on  a  resonator  duct  (Ward,  1997).  This  limitation 
makes  it  very  difficult  to  model  an  annular  prime  mover  and  incorporate  measured 
temperatures  along  the  duct  into  the  numerical  model. 

Because  of  the  disadvantages  of  applying  the  DeltaE  program  to  our  problem  a 
decision  was  made  to  develop  our  own  program.  This  program,  written  in  MATLAB,  was 
validated  by  comparing  the  results  for  a  simplified  annular  prime  mover  with  DeltaE 
results. 

3.    MATLAB  Program 

a.     Approach 

One  important  aspect  of  this  dissertation  is  to  develop  a  numerical  analysis 
program  that  is  more  tailored  to  the  annular  resonator  problem  than  DeltaE.  The  program 
implements  the  shooting  method  because  it  appears  to  be  a  straightforward  way  to  solve 
boundary  value  problems.  This  method  also  provides  a  systematic  approach  to  taking  a  set 
of  ranging  shots  that  allows  us  to  improve  our  aim  systematically  (Keller,  1992;  Press  et 
al.,  1992).  The  program  uses  functions  that  already  exist  in  MATLAB  where  possible. 
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A  complete  explanation  of  the  procedure  of  the  shooting  method  can  be  found  in 
Press  (1992)  and  Gerald  (1994).  Only  a  summary  is  presented  here.  Most  differential 
equations  of  order  higher  than  first  can  be  reduced  to  coupled  sets  of  first-order  differential 
equations.  The  two  point  boundary  value  problem  is  equivalent  to  obtaining  the  solution  to 
a  set  of  TV  coupled  first-order  differential  equations,  satisfying  N  periodic  boundary 
conditions  at  the  starting  point  xa  and  at  the  end  point  xb. 

We  start  by  writing  Eqs.  (2.25)  through  (2.29)  in  a  general  form 

^y^   =  /i(*  yv  Ja. >yN\  *  =  U N,  (2.30) 

ax 

with  the  periodic  boundary  conditions  at  the  starting  point  ;ca  and  final  point  jcb  expressed  as 

?*(*.»  yv  ^ >yN)  =  yi{xb>  yv  y^ >yN)>        i  =  12, ,  n,   (2.31) 

where,  in  our  problem,  N  equals  five. 

To  solve  this  boundary  value  problem,  an  initial  value  problem  is  created  by 
assuming  five  initial  values.  The  MATLAB  function  ode45,  which  implements  fourth  and 
fifth  order  Runge-Kutta  formulas  is  then  used  to  integrate  the  differential  equations  from  ;ca 
to  jcb.  The  local  error  term  for  the  fourth-order  Runge-Kutta  method  is  0(/z5)  and  the  global 
error  would  be  0(/i4),  where  h  is  the  step  size  of  the  integration  (Celia  and  Gary,  1992; 
Gerald  and  Wheatley  1994).  Now,  at  the  final  point  jcb,  a  discrepancy  vector  F  is  defined 
with  the  dimension  of  N,  whose  components  measure  how  far  the  first  solution  attempt  is 
from  satisfying  the  boundary  conditions  at  jcb.  The  problem  now  becomes  one  of  finding 
the  vector  F  such  that 

F,  (yv  y2 yN)  =  °>       i  =1-  2>->N-  (2-32) 
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To  find  the  roots  of  this  equation  the  Newton-Raphson  method  is  employed.  This  method 
provides  a  very  efficient  means  of  converging  to  a  solution,  if  a  sufficiently  good  initial 
guess  is  given  (Press,  1992). 

Letj>  denote  the  entire  vector  of  values  v..  In  the  neighborhood  of  y ,  each  of  the 
functions  Fy  can  be  expanded  in  a  Taylor  series  about  v0  as  (Press,  1992) 

FM  =  ?&<>+%)  =  FiM  +  Ix"  &7-  +0(6y2).  (2.33) 

The  matrix  of  partial  derivatives,  known  as  the  Jacobian  matrix,  J  is  defined  as 

7,  =   -A  (2.34) 

Since  it  is  difficult  to  compute  the  matrix  J  analytically  in  this  problem,  finite  differences 
are  used  to  compute  this  matrix.     In  the  MATLAB  program,  this  is  accomplished  by 
perturbing  each  yi  individually  and  finding  the  value  of  the  vector  of  AF-JAy-fc). 
Equation  (2.34)  can  be  expressed  in  the  matrix  form 

F(y  +  by)  =  F(y)  +  J  •by  +0(dy2)  .  (2.35) 

To  obtain  a  set  of  linear  equations  for  the  correction  vector  6y,  terms  of  order  by2 
and  higher  are  neglected.   Setting  F(y+$y)  =  0,  gives 

J.dy  =    -  F(v).  (2.36) 

Provided  J  is  nonsingular,  by  is  given  by 

fy  =   -J'F{y).  (2.37) 

In  practice,  J'1  is  obtained  by  LU  factorization.  In  MATLAB  notation,  by  is  given  by 
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by  =  -J\F(y),  (2.38) 

where  \  is  used  to  denote  left  division. 

The  LU  factorization  implements  Gaussian  elimination.  Because  it  involves  the 
least  amount  of  arithmetic  operations,  Gaussian  elimination  is  generally  considered  to  be 
one  of  the  most  efficient  computation  methods  in  the  problem  of  solving  a  system  of  linear 
equations  (Leon,  1994;  Lindfield  and  Penny,  1995).  It  is  worth  pointing  out  that  using 
J\F(y)  instead  of  inv(J)*F(y)  in  MATLAB  is  two  or  three  times  faster  and  produces 
residuals  on  the  order  of  machine  accuracy  relative  to  the  magnitude  of  the  data  (The  Math 
Works,  1994).  Once  by  is  found,  the  solution  vector  is  then  updated  as 

ynew  =  yM  +  &  (2-39) 

The  process  is  iterated  until  the  solution  converges,  i.e.  when  the  value  of  the  normalized 
norm  of  the  discrepancy  vector  F  is  smaller  than  a  desired  threshold. 

In  order  to  formulate  a  program  to  model  the  annular  prime  mover,  the  parameters 
to  be  updated  at  each  pass  of  the  integrations  must  be  identified.  These  parameters  are 
selected  to  be:  the  real  and  imaginary  parts  of  the  volume  velocity  U\  and  the  complex 
frequency  CO,  and  the  time-averaged  energy  flux  along  the  stack H2 .  Re(p,)  and  Im(p,)  are 
kept  constant  at  x  =  0,  which  simply  normalizes  the  amplitude  of  the  pressure  variation. 
In  fact,  for  given  values  of  Re(/?,)  and  Im(p,),  a  unique  solution  can  be  obtained  for  a 
certain  value  of  H2 .  The  periodic  boundary  conditions  are  satisfied  by  matching  the 
acoustic  pressure  p]  and  the  acoustic  volume  velocity  U\  between  the  starting  and  final 
points.  Also  the  temperature  at  the  hot  heat  exchanger  is  matched  to  a  certain  value  while 
keeping  the  temperature  of  the  ambient  heat  exchanger  a  constant. 
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b.      Description 

The  complete  MATLAB  program  is  described  in  Appendix  C.  This  section 
provides  a  brief  explanation  of  how  the  program  works.  First,  the  coordinate  system  for 
the  prime  mover  (as  is  shown  in  Fig.  2.2)  is  defined.  The  hot  heat  exchanger  is  at  one  end 
of  the  system.  The  end  of  the  hot  heat  exchanger  is  at  x  =  Leff.  The  choice  of  the  position 
of  x  =  0  is  arbitrary,  however  by  putting  the  stack  and  heat  exchangers  at  the  end,  only 
four  distinct  regions  are  required.  Had  the  stack  been  placed  in  the  middle,  there  would 
have  been  five. 


Duct 


Ambient  Heat 
Exchanger 


Prime  Mover    Hot  Heat 
Stack  Exchanger 


— f— 

Periodic 
Boundary 


3  Periodic 
j^  Boundary 


x=0 


x=Leff 


Figure  2.2    Coordinate  system  used  in  the  numerical  analysis  of  the  annular  prime 
mover. 


Initial  values  of  the  real  and  imaginary  parts  of  pv  a  measured  temperature  profile 
Tm(x)  along  the  duct  and  ambient  heat  exchanger,  and  a  temperature  for  the  hot  heat 
exchanger,  Th  are  given.  The  real  and  imaginary  parts  of  U\,  the  real  and  imaginary  parts 
of  the  resonance  frequency  co  and  the  time-averaged  energy  flux  along  the  stack  H2  are 
guessed.  It  is  assumed  that  the  temperature  of  the  duct  and  the  ambient  heat  exchanger  are 
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held  at  certain  values  by  external  means.  A  set  of  reasonable  initial  guesses,  which  can  be 
estimated  from  the  output  of  DeltaE,  is  the  key  to  successful  convergence  of  the  program. 
The  desired  temperature  of  the  hot  heat  exchanger  (Th  in  the  program)  is  accomplished  by 
adjusting  the  time-averaged  energy  flux  along  the  stack  H2 .  The  targets  to  be  matched  are 
the  real  and  imaginary  parts  of/?,  and  U\  at  x  =  0  and  x  =  Leff,  and  the  temperature  of  the 
hot  heat  exchanger.  This  matching  is  accomplished  by  adjusting  the  real  and  imaginary 
parts  of  Uu  the  real  and  imaginary  parts  the  complex  angular  frequency  0),  and  the  time- 
averaged  energy  flux  along  the  stack.  It  is  also  assumed  that  there  are  no  temperature 
gradients  in  the  ambient  and  hot  heat  exchanger. 

Once  the  required  input  parameters  are  determined.  The  program  solves  the 
system  of  equations  for  the  five  unknown  variables.  If  the  resulting  new  pv  U\,  and  Th 
differ  significantly  from  the  initial  guesses  it  is  necessary  to  update  the  guesses  and  do  the 
integration  again.  During  the  integration,  all  the  temperature  dependent  parameters  and  gas 
parameters  are  calculated  according  to  the  temperature  distribution.  This  whole  process  is 
repeated  to  convergence. 

c.     Validation 

This  program  was  validated  by  comparing  its  solution  to  that  using  DeltaE  for  a 
simple  case,  in  which  it  is  assumed  that  the  temperatures  of  the  entire  duct  and  the  ambient 
heat  exchanger  are  held  at  293  K  by  some  external  means.  The  desired  temperature  of  the 
hot  heat  exchanger  is  achieved  by  adjusting  the  quantity  H2 . 

The  comparisons  of  the  results  from  the  two  programs  are  shown  in  Figs.  2.3 
through  2.8.  Figure  2.3  shows  resonance  frequency  vs.  AT,  with  AT  increasing  from 
0  K  to  240  K.   Figure  2.4  shows  1/(2   vs.  AT  over  the  same  temperature  difference  span 
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as  in  Fig.  2.3.  The  overall  agreement  of  resonance  frequencies  and  l/Q  from  the  two 
methods  is  good.  In  fact,  the  normalized  mean  square  errors  of  the  resonance  frequency 
between  two  methods  are  1.48xl0"2%  and  6.10xl0"3%  for  the  low  and  high  frequency 
mode,  respectively.  The  normalized  mean  square  error  of  l/Q  between  two  methods  are 
0.684%  (low  frequency  mode)  and  0.854%(high  frequency  mode).  The  value  of  l/Q 
calculated  from  DeltaE  is  determined  from  a  pole-zero  analysis  of  the  frequency  response. 
The  value  of  Q  so  determined  is  sensitive  to  the  bandwidth  and  weighting  used  in  the  fit. 
The  error  bar  attached  to  the  AT  =  200  K  data  point  indicates  the  uncertainty  and  accounts 
for  most  of  the  disagreement  at  higher  values  of  AT.  One  reason  for  the  dependence  of  Q 
on  bandwidth  is  the  multimode  excitation.  A  driver  was  incorporated  in  the  DeltaE 
program  to  find  the  frequency  response,  from  which  we  obtain  the  required  modal 
information.  The  driver  locations  are  selected  such  that  at  AT  =  0  K  there  is  only  one 
mode  being  excited.  However,  as  AT  is  increased,  the  mode  shapes  can  be  altered  and  the 
previous  driver  position  may  excite  both  modes  simultaneously.  This  can  be  seen  in  Figs. 
2.5  and  2.6  which  show  the  frequency  response  of  the  low  mode  and  high  mode  of  the 
prime  mover  from  DeltaE  for  three  values  of  AT.  For  higher  the  AT,  the  effect  of  high 
mode  excitation  is  clearly  visible  as  a  small  bump  in  the  frequency  response  of  the  low 
mode.  However,  no  such  bump  is  observed  in  the  frequency  response  of  the  high  mode. 

Figures  2.7  and  2.8  show  the  mode  shapes  for  the  low  frequency  mode  with 
AT  =  0  K  and  AT  =  100  K.  Figures  2.9  and  2.10  show  the  mode  shapes  for  high 
frequency  mode  with  AT  =  0  K  and  AT  =  100  K.  In  these  figures,  the  solid  lines 
represent  the  results  of  the  MATLAB  program  and  the  open  circles  represent  the  results  of 
the  DeltaE  program.  The  mode  shapes  from  DeltaE  have  been  normalized  to  the  highest 
value  of  the  pressure  amplitude.  The  amplitude  of  the  mode  shape  from  the  MATLAB 
program  is  determined  from  a  least  squares  fit  to  the  DeltaE  results. 
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It  is  seen  that  one  of  the  advantages  of  the  MATLAB  program  is  the  smooth 
mode  shape  curve.  The  agreement  of  the  mode  shapes  are  good  in  general.  However  the 
agreement  of  the  low  frequency  mode  is  not  as  good  as  for  the  high  frequency  mode, 
particularly  at  high  temperature  difference.  The  reasons  for  this  discrepancy  at  higher 
values  of  AT  are  not  fully  understood.  It  should  be  pointed  out  that  the  two  programs  treat 
the  heat  exchangers  differently.  The  MATLAB  program  assumes  that  the  mean  gas 
temperature  at  the  heat  exchanger  equals  the  temperature  of  the  heat  exchanger.  However, 
DeltaE  imposes  a  thermal  impedance  so  that  the  gas  and  heat  exchanger  temperatures  differ. 

Figure  2.8  displays  one  interesting  feature.  Both  program  predict  a  low  standing 
wave  ratio.  This  feature  is  not  present  in  the  high  mode. 

Although  there  remain  small  disagreements  between  DeltaE  and  the  MATLAB 
program,  the  major  features  are  in  good  agreement.  In  the  chapters  to  follow,  the 
MATLAB  program  will  be  used  to  analyze  the  annular  prime  mover  under  more  realistic 
conditions.  The  major  difference  is  that  the  measured  temperature  profiles  are  incorporated 
in  the  program. 
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Calculated  resonance  frequencies  VS  DeltaT  for  the  annular  prime  mover 
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Figure  2.3    Comparison  of  the  calculated  resonance  frequencies  of  the  annular  prime 
mover  as  determined  from  the  MATLAB  program  and  DeltaE.   Symbols  are  explained 

in  the  legend.  Agreement  is  excellent  at  all  values  of  AT.  The  normalized  mean  square 
error  is  1.48xl0"2%  for  the  low  mode  and  6.10xl0"3%  for  the  high  mode. 
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Figure  2.4   Comparison  of  1/g  of  the  prime  mover  vs.  AT  as  determined  from  the 

MATLAB  and  DeltaE.  An  error  bar  is  indicated  for  the  low  mode  at  AT  =  200  K. 
Because  of  modal  overlap,  determination  of  \IQ  in  DeltaE  is  not  without  uncertainty. 
Agreement  is  better  at  low  AT.  The  normalized  mean  square  error  is  0.684%  for  the 
low  mode  and  0.854%  for  the  high  mode. 
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Frequency  response  of  the  low  mode  for  the  annular  prime  mover 
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Figure  2.5     Frequency  response  of  the  prime  mover  from  DeltaE  (low  mode)  at 

AT  =  0  K,     100  K    and  200  K.       Note   the  small  bump   at  about  436  Hz    on 
AT  =  100  K  and  200K  curves  which  is  a  result  of  mode  overlap. 
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Frequency  response  of  the  high  mode  for  the  annular  prime  mover 
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Figure  2.6     Frequency  response  of  the  prime  mover  from  DeltaE  (high  mode)  at 

AT  =  0  K,  100  K  and  200  K.   Notice  that  there  is  no  visible  bump  from  overlap. 
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Low  mode  with  DeltaT  =  0  Kelvin  ;  frequency  =  419.2+6.038i 
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Figure  2.7  Comparison  of  mode  shapes  of  the  prime  mover  at  AT  =  0  K    for  the  low 

frequency  mode  as  calculated  with  the  MATLAB  program  and  DeltaE.   The  beginning 
of  the  stack  region  is  indicated  by  the  dashed  vertical  line. 


40 


Low  mode  with  DeltaT  =  100  Kelvin  ;  frequency  =  422.3+6. 192i 
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Figure  2.8    Comparison  of  mode  shapes  of  the  prime  mover  at  AT  =  100  K    for  the 

low  frequency  mode  as  calculated  with  the  MATLAB  program  and  DeltaE.   Although 
both  programs  predict  the  same  general  mode  shape,  the  agreement  is  not  as  good  at 

AT  =  0  K. 
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High  mode  with  dT  =  0  Kelvin  ;  freq  =  437.1  +2.01 3i 


0.4 
Position  (m) 


Figure  2.9  Comparison  of  mode  shapes  of  the  prime  mover  at  AT  =  0  K  for  the  high 
frequency  mode  as  calculated  with  the  MATLAB  program  and  DeltaE. 
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High  mode  with  dT  =  100  Kelvin  ;  freq  =  437+2. 137i 


0.3  0.4  0.5 

Position  (m) 

Figure  2.10  Comparison  of  mode  shapes  of  the  prime  mover  at  AT  =  100  K    for  the 

high  frequency  mode  as  calculated  with  the  MATLAB  program  and  DeltaE.     The 

agreement  is  much  better  than  for  the  low  mode  at  AT  =  100  K. 
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E.    END  EFFECTS 

The  usual  boundary  conditions  applied  to  a  change  in  cross  section  are  continuity  of 
pressure  and  volume  velocity.  However,  at  the  discontinuity  an  additional  acoustic 
impedance  is  introduced  owing  to  the  abrupt  change  of  cross  sectional  area  (Miles,  1946 
and  1947;  Karal,  1953;  Pierce,  1991).  As  a  matter  of  fact,  we  have  shown  (Choe,  1997) 
that  it  is  critical  to  take  into  account  end  corrections  at  the  constriction  to  get  good 
agreement  between  the  measured  and  predicted  resonance  frequency  and  mode  shapes  for  a 
constricted  annular  resonator. 

Sudden  changes  in  duct  cross  section  produce  sudden  changes  in  acoustic  pressure 
amplitude.  Below  cutoff,  the  effect  is  equivalent  to  a  lumped-parameter  impedance.  For 
example,  a  constriction  in  a  duct  gives  rise  to  a  concentration  of  flow  through  the 
constriction  that  consequently  increases  the  kinetic  energy  density.  As  is  discussed  by 
Morse  and  Ingard  (1968),  if  the  net  increase  is  concentrated  in  a  region  small  compared  to  a 
wavelength,  the  total  increase  can  be  characterized  by  a  lumped  acoustic  inductance.  If  the 
constriction  also  produces  an  increase  in  viscous  energy  loss,  this  addition  can  be  modeled 
as  a  lumped  acoustic  resistance. 

This  equivalent  lumped  impedance  can  be  obtained  from  various  methods.  Only  two 
methods,  the  conformal  transformation  method  and  the  higher  order  mode  method,  are 
discussed  here. 

After  computing  the  analogous  acoustical  impedance  Z,  the  boundary  condition  for  the 
pressure  becomes. 

P2(D  =  px(D  -  Z*UX{1),  (2.40) 
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where  px  is  the  pressure  in  the  unconstricted  region  and  p2  is  the  pressure  in  the 
constriction.  Ux  is  the  volume  velocity  in  the  unconstricted  region  and  the  constriction 
junction  is  at  x  =  /.  The  volume  velocity  boundary  condition  is  unchanged. 

1.     Conformal  Transformation  Method 

The  first  method  used  to  compute  the  lumped  impedance  introduced  by  a  constriction 
is  discussed  by  Morse  and  Ingard  (1968).  To  obtain  the  analogous  impedance  of  a 
constriction,  they  compute  the  extra  kinetic  energy  and  viscous-energy  loss  produced  by 
steady  flow  through  the  constriction.  To  find  the  steady-state  flow,  the  interior  of  the  half- 
duct  is  first  transformed  to  the  upper  half  of  a  complex  w  plane  by  the  Schwartz-Christoffel 
transformation.  Next,  depending  on  the  geometry  of  a  specific  problem,  a  further 
transformation  is  required  to  go  from  the  w  plane  to  the  velocity-potential  plane  in  which 
the  lumped  impedance  can  be  readily  obtained. 

Morse  and  Ingard  applied  this  technique  to  the  case  of  a  rectangular  duct  of  depth  d 
and  width  b  for  x  <  0  and  a  for  x  >  0  (as  shown  in  Fig.  2.11).  The  lumped-circuit 
inductance  and  resistance  corresponding  to  the  sudden  increase  in  the  cross  section  of  a 
rectangular  duct  are 


L,  =     P- 


.2 


Ttd 
and 

p(o8v  a-b 


(a-b)       a  +  b  (a  +  b)' 

In +   In 


2ab         a  —  b  4ab 


(2.41) 


R,  = 


lad      b 


'         a2-b2    ,    a  +  b} 

1  +  In 

y  nab  a-b  j 


(2.42) 


Notice  that  both  L  and  R  vanish  when  a  =  b  (i.e.,  when  there  is  no  change  in  width). 
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Figure  2.11     Section  of  a  rectangular  duct  with  sudden  change  of  width  in  the  y 
direction. 


2.    Higher  Order  Mode  Method 

This  method  involves  including  higher  order  modes  near  the  constriction  ends 
(Muehleisen,  1996).  In  this  work,  Muehleisen  examines  the  reflection,  transmission, 
radiation,  and  coupling  of  higher  order  modes  at  a  discontinuity  in  finite  length  rigid  wall 
rectangular  ducts.  Using  a  method  of  generalized  scattering  coefficients,  an  analytic 
expression  for  the  reflection  and  transmission  of  higher  modes  at  size  discontinuities  is 
developed. 

From  the  scattering  matrix  formulation  an  expression  for  the  lumped  acoustic 
impedance  can  be  found  (Muehleisen,  1997).  This  method  does  not  include  losses, 
however,  and  so  the  resistance  obtained  by  conformal  mapping  is  added  to  the  impedance. 


46 


III.     EXPERIMENTAL  APPARATUS  AND  PROCEDURE 

A.  INTRODUCTION 

Descriptions  of  the  experimental  apparatus  and  data  acquisition  procedure  are 
presented  in  this  chapter.  The  chapter  begins  with  a  discussion  of  the  annular  resonator. 
This  resonator  was  used  in  preliminary  work  to  measure  the  modal  properties  of  a 
constricted  annular  resonator  in  the  absence  of  a  stack  (Choe,  1997).  This  discussion  is 
followed  by  a  detailed  description  of  the  annular  prime  mover.  Next,  the  constricted 
annular  prime  mover  is  discussed,  followed  by  a  description  of  how  the  prime  mover  is 
instrumented  with  microphones  and  thermocouples.  The  calibration  of  the  microphones  is 
discussed  in  this  subsection  as  well.  The  chapter  concludes  with  a  description  of  the 
procedure  followed  during  data  acquisition. 

B.  EXPERIMENTAL  APPARATUS 

1.    Annular  Resonator  System 

The  annular  resonator,  shown  in  Fig.  3.1  (without  the  top)  and  Fig.  3.2  (with  the  top) 
has  been  used  in  a  preliminary  work  to  investigate  the  modal  properties  of  a  constricted 
annular  resonator  (Choe,  1997).  The  inner  and  outer  radius  of  the  resonator  is 
10.00  ±  0.01  cm  and  15.30  ±  0.01  cm,  respectively.  The  outer  side  wall  is  a  1  cm  thick 
aluminum  ring  with  15.30  ±  0.01  cm  inner  radius.  A  solid  aluminum  10  cm  radius 
cylinder  that  is  concentric  with  the  outer  ring  forms  the  inner  side  wall.  The  height  of  the 
resonator  is  5.00  ±  0.01  cm.    The  dimensions  of  the  resonator  were  chosen  to  give  a 
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longitudinal  resonance  frequency  in  the  300  to  500  Hz  range,  which  is  well  below  the 
cutoff  frequency  for  the  first  cross  mode,  as  discussed  in  Chapter  II. 

Since  the  prime  mover  will  be  operated  at  temperatures  over  200  °C,  the  plexi-glass 
top  described  in  Choe  has  been  replaced  by  a  0.64  ±  0.01  cm  thick  aluminum  top  with  a 
radius  of  16.01  ±  0.01  cm.  The  top  is  bolted  with  a  10.7  cm  long  threaded  rod  at  the 
center  and  three  equally-spaced  7.1  cm  long  clamp  bars  on  the  edges.  Ten  microphones 
are  installed  in  the  top  to  measure  the  spatial  pressure  distribution  inside  the  driven 
resonator  below  onset  and  observe  the  behavior  of  the  sound  generated  above  onset.  To 
monitor  the  temperatures  in  the  resonator,  5  type-E  thermocouples  are  installed  through  the 
top  and  another  5  type-E  thermocouples  are  installed  in  the  stack/heat  exchanger  assembly. 
A  detailed  description  of  the  microphone  and  thermocouple  installation  will  be  presented 
later  in  this  chapter. 

It  is  very  difficult  to  move  the  stack  and  heat  exchanger  assembly  after  it  has  been  put 
into  the  resonator  because  of  the  complex  wiring  and  the  fragility  of  the  glass  components. 
Therefore,  three  holes  (3.50  mm  diameter)  were  drilled  in  the  bottom  of  the  resonator  to 
allow  the  driver  location  to  be  moved  relative  to  the  stack.  Depending  on  which  mode  is  to 
be  excited,  the  driver  was  then  connected  to  the  desired  hole  via  a  0.4  cm  O.D.,  3.5  cm 
long  plastic  tube.  The  two  unused  holes  were  plugged  with  thumbtacks  that  had  some 
Dow  Corning  high  vacuum  grease  applied  to  the  tips. 

To  ensure  tight  sealing  for  a  high  Q,  a  thin  layer  of  Dow  Coming  high  vacuum  grease 
was  applied  to  the  surface  of  the  aluminum  ring  and  the  center  piece  before  the  top  was  put 
on  and  bolted  down.  In  addition  to  the  bolt  in  the  center  and  the  three  clamps  mentioned 
earlier,  6  extra  equally  spaced  clamp  bars  are  applied  around  the  edge  of  the  top  to  ensure  a 
tight  seal. 
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The  instrumentation  is  shown  in  Fig.  3.3.  The  acoustic  signal  is  generated  by  a  40W, 
16  Q  external  compression  driver.  The  output  of  a  Hewlett  Packard  3562A  dynamic  signal 
analyzer  was  sent  through  a  Techron  5507  power  amplifier  to  the  driver.  The  signal 
analyzer  was  used  to  determine  the  frequency  response  of  the  resonator.  Typically,  the 
resonator  was  driven  with  a  random  noise  source  of  1  Vims  source  level.  The  pole-zero 
analysis  function  built  into  the  signal  analyzer  was  then  used  to  determine  the  Q  and  the 
resonance  frequency.  To  measure  the  mode  shape,  the  source  of  the  analyzer  was  changed 
to  a  lVrms  fixed  sinusoidal  signal  at  the  resonance  frequency  found  in  the  previous 
measurement. 

2.     Stack  Assembly 

The  prime  mover  stack  assembly  consists  of  the  stack  plates  and  the  stack  holder.  The 
first  challenge  was  to  select  a  suitable  material  for  stack.  The  stack  needed  to  withstand 
high  temperatures  and  the  differential  expansion  imposed  upon  it  from  the  high  gradients. 
Glass  has  a  low  thermal  conductivity  in  comparison  to  most  metals  (Johanson  Companies, 
1996)  and  is  able  to  withstand  higher  temperatures  than  most  plastics  and  other  polymers. 
An  AF-45  thin  borosilicate  glass  slide  (made  by  Precision  Glass  &  Optics)  was  chosen  to 
be  the  material  for  stack.  Some  properties  of  this  AF-45  glass  are  provided  in  Appendix  D. 
AF-45  is  a  modified  borosilicate  glass  with  a  high  content  of  BaO  and  A1203  and  features  a 
low  thermal  conductivity  and  a  low  coefficient  of  thermal  expansion  of  45.0  x  lO^K"1 
(20-300°  C).  The  glass  was  accurately  cut  to  a  size  of:  43.65+/-0.05  mm  wide  and 
23.90+/-0.10  mm  long  and  is  0.145+/-0.02  mm  thick. 

The  stack  holder  (detailed  dimensions  of  which  are  given  in  Appendix  E)  consists  of 
two  horizontal  bars  (5.28  cm  long)  connected  to  two  vertical  bars  (4.85  cm  long)  as  shown 
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in  Fig.  3.4.  Tongue-and-groove  joints  are  machined  in  the  ends  of  the  bars  to  ensure  that 
they  form  a  sturdy  frame.  Sixty-one  0.508  mm  deep,  0.203  mm  wide,  equally-spaced  slits 
were  cut  in  the  inside  of  each  of  the  two  horizontal  bars.  The  center-to-center  spacing 
between  the  slits  is  0.762  mm.  The  gaps  between  the  cover  glass  plates  are  0.610  mm. 
After  the  stack  holder  was  assembled,  61  glass  slides  were  then  carefully  inserted  into  the 
slits.  The  overall  porosity  for  the  prime  mover  stack  assembly,  including  the  frame,  is 
approximately  0.61. 

3.     Ambient  Heat  Exchanger  Assembly 

The  ambient  heat  exchanger  is  used  to  remove  heat  from  the  stack  and  to  maintain  one 
end  of  the  stack  close  to  ambient  temperature.  The  ambient  heat  exchanger  assembly  is 
constructed  as  shown  in  Fig.  3.5  and  the  detailed  dimensions  are  given  in  Appendix  E. 
The  heat  exchanger  is  designed  to  ensure  that  the  temperature  difference  between  the  center 
and  the  edge  of  the  heat  exchanger  is  less  than  one  degree  C  for  a  35  W  load.  The  ambient 
heat  exchanger  holder  was  made  in  the  same  way  as  the  stack  holder,  except  that  it  uses 
copper  instead  of  aluminum  and  has  only  10  slits.  Also,  a  central  bar  was  incorporated  to 
supplement  the  heat  transfer  between  the  heat  exchanger  and  the  resonator  body.  The  slits 
are  0.635  mm  deep,  0.71 1  mm  wide,  and  equally  spaced  at  4.064  mm  (center  to  center). 
After  the  holder  was  made,  twenty  Alloy-100,  2.30  mm  long,  5.10  mm  wide,  and  0.68 
mm  thick  copper  fins  were  cut  and  soft-soldered  into  the  slits  (10  on  either  side  of  the 
central  bar).  The  overall  porosity  for  the  ambient  heat  exchanger  is  0.62  (including  the 
frame).  The  next  step  was  to  soft-solder  a  5.256  cm  x  4.977  cm  rectangular  copper  woven 
wire  cloth  (20  x  20  mesh)  on  one  side  of  the  heat  exchanger  assembly.   The  function  of 
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this  wire  cloth  is  to  increase  the  thermal  contact  between  the  prime  mover  stack  and  the 
ambient  heat  exchanger. 

The  fin  spacing  is  much  larger  than  would  be  used  in  a  practical  prime  mover.  The 
acoustical  properties  of  an  constricted  annular  resonator  are  very  sensitive  to  both  the 
porosity  and  flow  impedance  of  the  constriction.  It  was  decided  to  make  the  flow 
impedance  of  the  heat  exchangers  small  compared  to  that  of  the  stack.  This  choice 
necessitated  using  large  fin  spacing,  which  reduced  the  effectiveness  of  the  heat  exchanger. 

4.    Hot  Heat  Exchanger  Assembly 

The  function  of  the  hot  heat  exchanger  is  to  supply  heat  to  the  stack.  It  consists  of 
nickel-chromium  wire  woven  on  a  non-conducting  frame.  Since  nickel-chromium  wire  is 
used  to  provide  heat  to  the  system,  the  hot  heat  exchanger  has  to  provide  electrical 
insulation  to  the  aluminum  resonator.  Because  of  the  need  for  electrical  insulation  and  high 
temperature  operation,  virgin  electrical  grade  PTFE  (Teflon)  is  chosen  for  the  hot  heat 
exchanger  frame.  The  hot  heat  exchanger  assembly  is  shown  in  Fig.  3.6.  Drawings  of 
fabricated  parts  are  provided  in  Appendix  E.  The  Teflon  frame  has  outer  dimensions  of 
4.978  cm  wide,  5.283  cm  high,  and  0.813  cm  thick.  The  inside  of  this  frame  was  milled 
out  to  a  size  of  4.470  cm  wide  and  4.663  cm  high.  Seventeen  holes  are  drilled  onto  the 
two  vertical  sides  with  an  equal  spacing  of  0.521  cm  (8  on  one  side  and  9  on  the  other 
side).  Into  each  of  the  holes  is  inserted  a  0-80  stainless  steel  screw.  A  92.1  cm,  #  28  gage 
nickel-chromium  wire  with  a  total  resistance  of  12  Q.  was  wound  around  the  17,  0-80 
screws.  To  avoid  direct  contact  between  the  nickel-chromium  wire  and  the  stack,  stainless 
steel  flat  washers  are  used  on  the  screws  for  spacers.  The  ends  of  the  nickel-chromium 
wire  are  wrapped  with  heat  shrinkable  flexible  PVC  tubing  to  insulate  the  wire  from  the 
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resonator.  The  wire  is  then  fed  through  Teflon  tubing  which  is  placed  and  epoxied  in  the 
bottom  wall  of  the  resonator. 
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Figure  3.1  The  annular  resonator  (without  the  top). 
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Figure  3.2  The  annular  resonator  (with  the  top). 


54 


S 

s 

>H 

o 

-J 

UJ 

-j 

0- 

D 

00 

£ 

Oi 

uj 

CO 
on 

O 

(N 

a. 

VD 

Cu 

sc 

I 


oi 

00   ^ 

00 
UJ 

a. 

D 

£  o 

O 

£  £ 

u.. 

o  oi 

0 

oi 

fe    ^ 

tq 

So 
I  z 

£  z 

uj  z 

2  < 

o 

o 

en 

B" 

Oi 

n 

00 

o 

Qu 

UJ 

p 

z 

b 

uj 

*\   8  p 

o 

X 

j 

b  ^■^v 

0. 

UJ 

/  ^      \ 

o  . 

00           y 

UJ 

->  J  ®    J 

oi 

u 

z 

o 

a. 

o 

2 
o 

oi 
U 

p" 

u 

s 

b 

ft 

o 

1 


u  oi 

UJ 


m — 
y 

Z  H 

el 


cT 

i  1 

tu 

in 

D- 

y 

o 

z 

< 
DC 

u 

00 

O 

U 

_J 

o 

00 

O 

VO 

en 

m 

i — i 

CN 

_l 

X 

UJ 

y 

z 

z 

o 

< 

Oi 

u 

P 

Figure  3.3  Block  Diagram  of  Instrumentation  Setup. 
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Figure  3.4  The  prime  mover  stack  assembly. 
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Figure  3.5  The  ambient  heat  exchanger  assembly. 
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Figure  3.6  The  hot  heat  exchanger  assembly. 
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5.     Stack/Heat  Exchanger  Assembly 

The  stack  assembly,  the  ambient  heat  exchanger  assembly,  and  the  hot  heat  exchanger 
assembly  were  connected  with  four  5.08  cm  long,  00-90  stainless  steel  threaded  rods  on 
the  two  horizontal  sides  (as  shown  in  Fig.  3.7).  Teflon  spacers  were  used  to  ensure  proper 
spacing  between  the  stack  and  the  heat  exchangers.  The  spacing  between  the  glass  plate 
and  the  nickel-chromium  wire  is  approximate  1  mm.  Before  the  stack/heat  exchanger 
assembly  was  placed  into  the  resonator,  a  piece  of  thin  Teflon  tape  was  wrapped  on  the 
side  of  the  hot  heat  exchanger  to  avoid  electrical  short  circuit.  A  thin  layer  of  silicone  heat 
sink  compound  was  applied  to  the  side  of  the  ambient  heat  exchanger  frame  to  supplement 
the  heat  transfer  between  the  ambient  heat  exchanger  and  the  resonator  wall  and  to  seal  any 
gaps  between  the  frame  and  the  resonator. 
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Figure  3.7  The  stack/heat  exchanger  assembly. 
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6.     Constriction 

Placing  another  constriction  in  the  resonator  provides  the  possibility  for  an  annular 
prime  mover  to  reach  onset.  Such  constrictions  were  previously  used  in  related  work  by 
Choe  (Choe,  1997).  The  constrictions  are  made  of  PVC  (Ploy vinylchloride)  because  of  its 
rigidity  and  ease  of  fabrication.  The  area  ratio  is  defined  by  the  ratio  of  open  area  in  the 
constriction  (which  is  the  shaded  block  in  Fig.  3.8)  to  the  cross-sectional  area  of  the 
resonator.  In  Fig.  3.8,  S^  stands  for  the  open  area  and  S2  stands  for  the  open  cross- 
sectional  area  of  the  constriction.  The  table  in  Figure  3.8  gives  the  geometrical  parameters 
of  the  constrictions.  Constrictions  were  used  having  area  ratios  of  0.7,  0.3  and  0.1  with  a 
total  length  of  45°.  The  45°  length  was  achieved  by  stacking  two  1 1.25°  constrictions  and 
one  22.5°  constriction  together.  Vacuum  grease  was  applied  to  the  joining  sides  of  the 
individual  constrictions  to  ensure  a  tight  fit. 
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Figure  3.8(a)  Nominal  geometry  of  the  constrictions. 
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Figure  3.8(b)  Measured  geometry  of  the  constrictions. 

7.     Microphones 

In  this  experiment,  8  Panasonic  electret  WM-60AT  microphones  and  two  ENDEVCO 
piezoresistive  pressure  transducers  (Model  #8510B-5,  Serial  #  G63T,  sensitivity  53.85 
mV/psi,  and  Model  #8530C-15,  Serial  #  AGLC2,  sensitivity  10.89  mV/psi)  are  used  to 
measure  the  pressure  distribution  inside  the  resonator.  The  two  ENDEVCO  piezoresistive 
pressure  transducers  are  intended  to  be  used  to  observe  the  behavior  of  the  sound  generated 
above  onset,  in  case  the  Pansonic  microphones  saturate.  The  Panasonic  microphones  are  6 

mm  in  height  and  5  mm  in  diameter,  with  a  sensitivity  of  -44  dB  ±  3  dB  re  1  V/Pa.    A 

Hewlett  Packard  6237  Triple  Output  Power  Supply  provided  9V  DC  to  each  Panasonic 
microphone  via  a  custom  made  microphone  selector. 

As  shown  in  Fig.  3.9,  the  8  Panasonic  microphones  were  spaced  by  45  degrees  in 
pre-drilled  holes  (0.605  cm  in  diameter,  0.50  cm  deep).  They  are  glued  into  the  holes  with 
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silicone  adhesive.  At  the  bottoms  of  these  holes,  0.82  mm  diameter  holes  were  drilled  to 
expose  the  microphones  to  the  resonator  cavity.  The  two  ENDEVCO  piezoresistive 
pressure  transducers  were  flush  mounted  in  the  cap  through  spacers  with  0.404  cm  in 
diameter  threaded  holes  in  the  center. 


0  degrees        45  degrees 


mic  #3 


mic  #1 


mic  #2  /     # 
270  degrees 


ENDEVCO  Pressure 
Transducers 


PANASONIC 
Microphones 

mic  #7 

90  degrees 


mic  #6 


180  degrees 


Figure  3.9  Installation  of  10  microphones  on  the  resonator  top. 


Prior  to  the  experiment,  the  8  Panasonic  microphones  were  calibrated  relative  to 
Microphone  #  1.  The  calibration  apparatus  is  shown  in  Fig.  3.10.  A  calibrating  chamber 
was  made  from  a  19.6  mm  long,  21.55  mm  ID,  and  27.22  mm  OD  circular  PVC  pipe.  A 
Panasonic  micro-speaker  (Model  #EAS-2P106C)  was  epoxied  into  one  end  of  the  pipe. 
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The  resonator  top  was  removed  and  the  calibration  chamber  was  positioned  over  each 
microphone.  A  thin  layer  of  Dow  Corning  high  vacuum  grease  is  applied  on  the  edge  of 
the  chamber  to  ensure  a  tight  sealing.  The  frequency  responses  were  then  measured  for 
each  microphone  over  a  300  ~  500  Hz  frequency  range  using  a  Hewlett  Packard  35665A 
Dynamic  Signal  Analyzer.  For  the  calibration,  the  signal  analyzer  was  set  to  swept  sine 
mode  with  a  integration  time  of  20  cycles,  a  settling  time  of  20  cycles,  and  a  resolution  of 
200  pts/sweep. 
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CHANNEL  2 
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b  5   6~ 
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Figure  3.10  Calibration  of  the  Panasonic  microphones  using  a  small  chamber. 

Microphone  #1  was  selected  as  the  reference  for  the  comparison  calibration.  The 
frequency  response  of  microphone  #1  when  covered  by  the  calibration  chamber  is  shown 
in  Fig.  3.1 1.  Each  microphone's  frequency  response  was  divided  by  that  of  microphone 
#1  to  create  a  correction  factor.   A  typical  curve  for  the  correction  factor  is  shown  in  Fig. 
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3.12.  After  finishing  the  calibrations,  one  more  measurement  was  made  on  Microphone  #1 
to  ensure  the  reproducibility  of  the  data.     With  this  particular  calibration  chamber,  a 

difference  of  less  than  1%  was  achieved.   The  correction  factor  for  the  8  microphones  is 

usually  in  the  range  of  0.75  to  1.75. 

When  measuring  the  mode  shape,  the  output  from  each  microphone  was  divided  by 
the  appropriate  correction  factor  to  obtain  the  corrected  amplitude  at  a  particular  location  and 
driving  frequency.  A  custom  made  microphone  switch  board  selects  each  microphone  to 
measure  the  pressure  amplitude  at  different  locations.  The  microphone  output  amplitude  is 
measured  using  a  Standford  Research  530  Lock-in  amplifier. 
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Figure  3.11     Frequency  response  of  microphone  #1   when  covered  by  the  small 
calibration  chamber. 
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Figure  3.12  Plot  of  correction  factor  of  microphone  #2  referenced  to  microphone  #1. 
The  correction  factor  at  400  Hz  is  0.940. 


8.     Thermocouple  Instrumentation 

Ten  Omega  type  E  (Model  #  5TC-GC-E-30-72)  chromel-constantan  0.254  mm 
diameter  glass  braid  insulated  thermocouples  are  used  to  monitor  the  temperatures  of  the 
stack,  the  heat  exchanger,  and  inside  the  resonator.   Type  E  thermocouples  were  selected 

over  other  types  because  they  provide  the  greatest  sensitivity  (in  V/°C)  and  because  of  their 

useful  temperature  range. 

As  shown  in  Fig.  3.13,  five  thermocouples  are  used  to  monitor  the  temperature  of  the 
stack  and  the  ambient  heat  exchanger.  Two  thermocouples  each  are  placed  on  the  hot  and 
ambient  sides  of  the  center  glass  slide,  one  in  the  center  and  the  other  at  the  edge.  These 
thermocouples  are  fed  through  1 .4  mm  diameter  holes  on  the  copper  woven  wire  cloth  of 
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the  ambient  heat  exchanger.  Another  thermocouple  (thermocouple  #  6  in  Fig.  3.13)  was 
placed  on  one  of  the  fins  of  the  ambient  heat  exchanger.  The  thermocouples  were  attached 
to  the  glass  plate  using  JB  weld  adhesive  because  its  high  temperature  capability.  The  size 
of  the  glue  spots  is  approximate  0.01  mm2.  These  glue  spots  are  approximate  0.5  mm 
from  the  edge  the  glass  plate.  The  wires  of  the  five  thermocouples  were  individually 
threaded  through  1 .4  mm  diameter  holes  in  the  resonator  bottom  adjacent  to  the  ambient 
heat  exchanger.  The  other  five  thermocouples  were  used  to  sense  temperature  inside  the 
resonator,  as  is  shown  in  Fig.  3.14.  They  were  fed  through  1.4  mm  diameter  holes  in  the 
top  and  positioned  approximately  3  cm  inside  the  resonator.  All  the  holes  in  the  resonator 
were  then  filled  with  BIPAX  TRA-BOND  epoxy  to  ensure  a  tight  seal. 

The  temperature  difference  across  the  stack  is  adjusted  by  the  voltage  supplied  to  the 
hot  heat  exchanger  which  is  controlled  by  a  Power  Design  3650R  DC  power  supply.  The 
temperatures  inside  the  resonator  and  in  the  stack  are  measured  by  a  Keithley  740  system 
scanning  thermometer. 
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Figure  3.13  Thermocouple  placement  on  the  stack  and  the  ambient  heat  exchanger. 
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Figure  3.14  Thermocouple  placement  in  the  resonator  top. 
C.    EXPERIMENTAL  PROCEDURE 

1.     Resonator  Setup 

The  stack/heat  exchanger  assembly  is  placed  at  a  fixed  location  relative  to  the 
microphones  and  thermocouples.  For  the  experiment,  three  driver  positions  are  used:  45°, 
90°,  and  180°  from  the  center  line  of  the  stack/heat  exchanger  assembly.  Unused  driver 
holes  are  plugged  with  thumbtacks  with  some  vacuum  grease  applied  to  the  tips.  The 
resonator  top  is  oriented  so  that  0°  is  set  at  the  stack  assembly  center  line.  The  0°  setting  is 
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performed  by  rotation  of  the  top  by  hand.  To  determine  the  error  in  positioning  the  top,  the 
top  was  positioned  five  times.  It  was  found  that  the  error  by  rotation  was  less  than  ±  1 .0°. 
So,  with  the  mechanical  error  ±  0.25°  in  drilling  the  microphone  holes,  the  total  error  of 
microphone  positioning  is  less  than  ±  1.25°.  With  the  stack/heat  exchanger  assembly 
positioned,  the  proper  amount  of  vacuum  grease  is  applied  on  the  open  surface  of  the 
resonator  and  the  top  is  secured  with  the  clamp  bars. 

For  the  measurements  with  a  constriction  in  the  prime  mover,  a  45°  long  constriction 
is  inserted  at  a  location  that  is  90°  from  the  center  line  of  the  stack/heat  exchanger  assembly 
to  the  center  of  constriction. 

2.     Determining  Resonator  Characteristics 

After  the  resonator  was  set,  the  Hewlett  Packard  3562A  Dynamic  Signal  Analyzer  is 
used  to  characterize  the  resonator.  The  swept  sine  mode  of  the  analyzer  provides  a  sine 
wave  signal  to  the  system  at  a  given  range  of  frequencies.  The  frequency  response  of  the 
system  is  used  to  determine  the  resonance  frequency  and  Q  of  several  modes  of  the 
resonator.  When  only  one  mode  is  excited,  the  peak  pressure  amplitude,  frequency,  and 
quality  factor  are  same  as  the  those  given  by  pole-zero  analysis.  When  two  modes  are 
excited  simultaneously,  pole-zero  analysis  is  used  to  extract  the  resonance  frequency  and 
quality  factor.  The  signal  analyzer  gives  poles  in  the  form  a+jb.  The  resonance  is  simply 
the  value  of  b  and  the  quality  factor  can  be  obtained  by  b/2a.  This  method  has  been 
explained  by  Choe  (1997).  These  measurements  are  performed  for  (1)  annular  prime 
mover  without  a  constriction,  and  (2)  annular  prime  mover  with  three  different 
constrictions  (45°  long,  porosities  of  0.7,  0.3  and  0.1)  placed  at  90°  relative  to  the  center 
of  the  stack/heat  exchanger  assembly.  The  driver  in  the  first  experiment  was  located  at  90° 
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(to  excite  the  low  frequency  mode)  or  180°  (to  excite  the  high  frequency  mode)  from  the 
center  of  the  stack  assembly.  In  experiment  (2),  the  driver  was  placed  at  45°  from  the 
center  of  the  stack  assembly,  between  the  stack  and  constriction. 

3.    Data  Collection 

Data  collection  is  performed  in  sets  of  annular  prime  mover  without  and  with  a 
constriction.  Before  data  collection,  the  output  of  the  Hewlett  Packard  6237B  Triple  output 
power  supply  is  set  to  DC  9V.  The  output  voltage  of  the  Power  Design  3650R  power 
supply  is  also  adjusted  to  a  certain  value  and  the  temperature  at  the  hot  end  of  the  stack  is 
allowed  to  reach  steady  state  (temperature  variation  <  0.1°C/10  sec).  The  measurement  is 
commenced  by  recording  the  temperatures  measured  by  the  Keithley  740  System  Scanning 
thermometer.  The  resonance  frequency  and  quality  factor  is  measured  by  the  Hewlett 
Packward  3562A  Dynamic  Signal  Analyzer  and  the  microphone  switch  board  is  set  to  the 
microphone  that  provides  the  maximum  output.  The  procedure  used  to  acquire  the 
resonance  frequency  and  quality  factor  was  described  in  subsection  1 .  After  finding  the 
resonance  frequency,  a  data  run  is  commenced  by  changing  the  source  of  the  Hewlett 
Packward  3562A  Dynamic  Signal  Analyzer  to  a  1  Vpp  fixed  sine  with  a  driving  frequency 
equal  to  the  resonance  frequency  obtained  by  previous  procedure.  The  acoustic  amplitude 
is  measured  with  each  microphone.  With  the  driver  running,  the  switch  board  selects  a 
particular  microphone  and  the  value  of  amplitude  at  each  position  is  measured  with  the 
Stanford  SR530  Lock-In  amplifier.  The  Lock-In  amplifier  is  set  to  magnitude-phase 
display  mode  with  a  time  constant  of  1  second.  At  this  time,  the  temperatures  are  recorded 
again  and  same  measurement  are  repeated  again.  Before  the  entire  set  of  experiments  is 
finished,  the  final  temperatures  are  recorded.  Then,  the  output  voltage  of  the  Power  Design 
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3650R  power  supply  is  changed  to  adjust  the  hot  end  temperature  and  the  same  procedure 
is  followed.  All  data  recorded  from  the  8  microphones  are  scaled  by  the  correction  factor 
obtained  from  the  microphone  calibration.  The  corrected  data  are  then  entered  in  a 
computer  program  to  plot  the  mode  shape.  This  concludes  data  acquisition  of  one  mode. 
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IV.     RESULTS  AND  DISCUSSION 

This  chapter  shows  the  comparison  between  theoretical  and  measured  results  for  two 
series  of  experiments  on  the  annular  prime  mover.  First,  measured  temperatures  in  the 
prime  mover  are  presented.  Next,  results  are  presented  for  resonance  frequency,  \IQ  vs. 
AT,  and  mode  shape  for  the  unconstricted  annular  prime  mover.  Results  for  the  annular 
prime  mover  with  a  constriction  located  90°  from  the  center  of  the  stack/heat  exchanger 
assembly  are  then  discussed.  Finally,  an  analysis  of  the  potential  for  onset  in  a  two-stack 
annular  prime  mover  is  presented. 


A.    ANNULAR  PRIME  MOVER 

The  measured  resonance  frequencies  and  Q's  are  determined  as  functions  of  AT  from 
pole-zero  analysis  of  the  frequency  response.  The  theoretical  resonance  frequencies  and 
quality  factors  are  determined  from  the  complex  frequency  calculated  using  the  measured 
temperature  distribution  in  the  prime  mover  as  an  input.  Figures  4.1  and  4.2  show  the 
measured  temperatures  as  functions  of  the  heater  voltage  at  four  locations  on  the  center 
plate  of  the  stack  for  the  low  mode  and  high  mode  of  the  prime  mover. 

Figure  4.1  shows  that  the  transverse  (center-to-edge)  temperature  difference  on  the 
ambient  side  of  the  stack  is  at  most  5.4  K.  However,  there  is  a  larger  transverse 
temperature  profile  on  the  hot  side.  At  the  highest  heater  voltage  the  hot-side,  center-to- 
edge,  temperature  difference  is  14.8  K.  The  transverse  temperature  profile  leads  to  some 
uncertainty  in  the  value  used  for  AT.  The  effects  of  a  transverse  AT  profile  on  the 
performance  of  a  stack  is  difficult  to  access.  Therefore,  the  transverse  profiles  will  be  used 
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to  place  error  bars  on  the  AT  values.  Figure  4.2  shows  the  center-to-edge  temperature 
difference  for  the  high  mode.  It  is  seen  that  the  ambient-side  temperatures  are  essentially 
the  same  as  in  Fig.  4.1.  However  the  hot  side  temperatures  show  a  larger  center-to-edge 
difference  than  in  Fig.  4. 1 .  The  difference  is  attributed  to  the  fact  that  the  resonator  top  had 
to  be  removed  to  reposition  the  driver  to  switch  from  exciting  the  low  mode  to  the  high 
mode.  It  is  entirely  possible,  and  observed,  that  while  removing  and  replacing  the  top  the 
stack/heat  exchanger  assembly  was  disturbed.  This  could  have  slightly  altered  the  hot  heat 
exchanger/stack  spacing,  resulting  in  the  observed  differences. 

In  the  results  to  follow,  all  the  low  mode  data  were  taken,  then  the  driver  was 
repositioned  and  all  the  high  mode  data  taken.  Therefore,  the  temperatures  displayed  in 
Figs.  4.1  and  4.2  should  be  representative.  Typical  measured  temperature  distributions  in 
the  prime  mover  are  shown  in  Figs.  4.3  (the  low  mode)  and  4.4  (the  high  mode).  The 
measured  temperature  profile  in  the  duct  (0  to  0.792  m)  and  the  ambient  heat  exchanger 
temperature  are  used  as  inputs  into  the  MATLAB  program.  The  measured  hot  end  stack 
temperature  is  a  target. 

Figure  4.5  shows  the  comparison  of  the  measured  and  computed  resonance  frequency 
for  the  low  mode  and  high  mode.  The  driver  location  is  90°  from  the  center  of  the 
stack/heat  exchanger  assembly  for  the  low  frequency  mode  and  180°  for  the  high  frequency 
mode.  It  is  seen  that  there  is  good  agreement  between  the  measured  and  predicted 
resonance  frequencies.  The  biggest  differences  between  the  measured  and  predicted 
resonance  frequencies  are  approximately  0.86  %  and  1.77  %  for  the  low  mode  and  high 
mode,  respectively. 

Figure  4.6  shows  \IQ  of  the  low  and  high  mode  for  the  annular  prime  mover  verses 
the  temperature  differences  across  the  stack.  In  general  the  agreement  is  good,  although 
there  is  a  tendency  to  under  predict  the  \IQ  (or  over  predict  Q),  indicating  the  presence  of 
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unaccounted  losses.  For  the  low  mode,  \IQ  decreases  as  the  temperature  difference  is 
increased  while  for  the  high  mode  \IQ  increases  as  the  temperature  increases.  It  is  evident 
from  these  results  that  the  annular  prime  mover  will  not  reach  onset  under  readily  attainable 
conditions.  As  a  matter  of  fact,  an  extrapolation  of  the  \IQ  curve  for  the  low  mode  gives  an 
onset  temperature  difference  of  approximately  1400  K. 

The  larger  discrepancies  of  the  resonance  frequency  and  UQ  for  the  high  mode  at  high 
temperature  gradient  can  be  caused  by  a  misjudging  of  the  stack  temperature.  In  the 
original  calculations,  the  temperatures  of  the  hot  (Th)  and  ambient  (Ta)  sides  of  the  prime 
mover  stack  are  calculated  by  a  weighted  average  of  the  temperatures  at  the  edge  and  center 
of  either  side  of  the  stack 


2T     +  T 
Th  =    Zite  +  ite  (4.1) 


and 


2T    +  T 

Ta  =       *         ae  ,  (4.2) 


where  Thc  and  rhe  are  the  center  and  edge  temperature  of  the  hot  side  of  the  stack, 
respectively.  Tac  and  Tat  are  the  center  and  edge  temperature  of  the  ambient  side  of  the  stack, 
respectively.  If  the  simple  average  of  Th  and  Ta  are  used  in  the  computations  as  opposed  to 
the  weighted  average,  the  resonance  frequencies  are  almost  unchanged  (difference  < 
0.03%).  This  negligible  change  is  to  be  expected  since  the  stack/heat  exchanger  assembly 
is  very  short  compared  to  the  effective  length  of  the  prime  mover 
(L  stack  assembly  I  ^  eff  =  0038).  As  a  results,  the  vertical  error  bars  in  Fig  4.5  are  smaller 
than  the  symbols  used  to  plot  the  data.  However,  l/Q  is  more  sensitive  to  uncertainties  in 
Th  and  Ta  than  is  the  resonance  frequency  since  all  thermoviscous  effects  are  most 
important  in  the  stack  region.   The  sensitivity  to  uncertainties  in  AT  is  represented  by  the 
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vertical  error  bars  shown  in  Fig.  4.6.  There  is  also  an  uncertainty  in  the  calculation  of 
temperature  difference  across  the  stack.  As  shown  in  Fig.  4.1  and  4.2,  the  temperatures  at 
the  center  of  the  hot  end  and  ambient  end  differ  from  those  at  the  edge.  The  largest 
differences  at  the  hot  end  are  14.8  K  and  33.3  K  for  the  low  and  high  mode,  respectively 
(heater  at  23  V  in  both  cases).  The  differences  of  the  center  and  edge  temperatures  in  the 
ambient  side  are  approximately  6  K  for  both  cases.  This  spread  is  represented  by 
horizontal  error  bars  shown  in  Figs.  4.5  and  4.6.  The  other  possible  source  of  the 
discrepancies  in  resonance  frequency  and  \IQ  at  high  AT  for  the  high  mode  is  the  entire 
temperature  distribution  of  the  prime  mover.  As  shown  in  Figs.  4.3  and  4.4,  the  major 
change  in  temperature  profile  is  concentrated  at  the  region  adjacent  to  the  hot  heat  exchanger 
and  stack.  It  is  difficult  to  accurately  characterize  the  entire  temperature  distribution  of  the 
prime  mover  using  only  10  thermocouples.  This  is  an  important  difference  between  the 
conventional  and  annular  prime  movers.  In  a  conventional  prime  mover,  there  is  typically  a 
hot  temperature  duct  and  an  ambient  temperature  duct.  In  an  annular  prime  mover  the  two 
temperatures  must  merge  along  the  duct. 

Figures  4.7(a),  (b),  and  (c)  show  the  mode  shapes  for  the  low  mode  at  different 
temperature  differences.  The  stack/heat  exchanger  assembly  is  located  between 
x  =  0.768  m  and  x  =  0.792  m,  indicated  by  the  dashed  line.  In  this  experiment,  the 
driver  is  located  90°  from  the  center  of  the  stack  which  corresponds  to  a  position  of 
0.183  m.  The  driving  frequencies  are  the  resonance  frequencies  as  computed  by  the 
MATLAB  program  and  the  AT's  are  calculated  using  Eqs.  4. 1  and  4.2. 

Figure  4.7(a)  shows  the  mode  shape  for  the  low  mode  for  AT  =  0  K.  In  this  case,  the 
stack/heat  exchanger  assembly  acts  simply  like  a  constriction  and  there  is  a  velocity 
antinode  near  the  center  of  the  stack/heat  exchanger  assembly.  The  agreement  between  the 
measured  and  calculated  mode  shape  is  good. 
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Figures  4.7(b)  and  (c)  show  the  mode  shapes  for  the  low  mode  at  temperature 
differences  of  80  K  and  200  K,  respectively.  The  discrepancy  between  the  measured  and 
computed  mode  shape  increases  as  the  temperature  difference  increases.  The  discrepancy 
appears  as  a  shift  in  mode  shape.  However,  in  both  cases,  the  measured  and  predicted 
mode  shapes  show  the  same  general  features.  In  particular,  the  standing  wave  ratio 
decreases  with  increasing  AT.  Even  with  this  obvious  diminishing  in  the  standing  wave 
ratio,  \IQ  only  changes  by  15%  from  AT  =  0  K  to  AT  =  220  K. 

Figure  4.8(a)  show  the  mode  shape  for  the  high  mode  with  AT  =  0  K.  The  agreement 
between  the  measured  and  calculated  mode  shape  is  excellent.  As  one  might  expect,  there 
is  a  pressure  antinode  near  the  center  of  the  stack/heat  exchanger  assembly.  Figures  4.8(b) 
and  (c)  show  the  mode  shapes  for  high  mode  at  temperature  difference  of  76  K  and  192  K, 
respectively.  The  agreement  between  measured  and  calculated  mode  shapes  are  better  than 
those  of  the  low  mode.  In  contrast  to  the  low  mode  shapes,  the  high  mode  shapes  do  not 
change  appreciably  as  AT  increases.  In  particular,  the  standing  wave  ratio  remains  high. 
This  difference  between  the  two  modes  may  be  caused  by  the  discontinuity  in  the  acoustic 
impedance  near  the  end  of  the  stack.  For  the  low  mode,  there  is  a  velocity  antinode  near 
the  center  of  the  stack  which  results  in  a  low  acoustic  impedance.  However,  the  acoustic 
impedance  of  the  stack  is  much  larger  for  the  high  mode.  As  a  result,  as  AT  increases,  the 
impedance  discontinuity  should  have  more  effect  on  the  low  mode  than  the  high  mode. 

The  most  important  result  of  this  section  is  that  the  annular  prime  mover  will  not  reach 
onset  under  easily  attainable  conditions  The  reason  is  that  neither  mode  shape  supports 
thermoacoustic  growth.  In  other  words,  the  low  mode  has  a  velocity  antinode  centered  on 
the  stack/heat  exchanger  assembly  whereas  the  high  mode  has  a  pressure  antinode  centered 
on  the  stack/heat  exchanger  assembly.  It  was  also  found  that  the  standing  wave  ratio  of  the 
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low  mode  diminishes  as  AT  increases.   However,  this  alteration  in  mode  shape  does  not 
correlate  with  increased  attenuation  in  the  prime  mover. 
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Temperature  of  the  prime  mover  stack  for  low  mode 
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Figure  4. 1    Temperature  of  the  prime  mover  stack  vs.  the  heater  voltage  (low 
mode). 
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Temperature  of  the  prime  mover  stack  for  high  mode 
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Figure  4.2    Temperature  of  the  prime  mover  stack  vs.  the  heater  voltage  (high 
mode). 
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Temperature  distribution  of  the  low  mode  of  the  annular  prime  mover 
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Figure  4.3    The  measured  temperature    distribution  of  the  prime  mover  (low 
mode,  heater  at  23  V).  The  stack  is  located  to  the  right  of  the  dotted  line. 


81 


Temperature  distribution  of  the  high  mode  the  annular  prime  mover 
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Figure  4.4    The  measured  temperature  distribution  of  the  prime  mover  (high 
mode,  heater  at  23  V).  The  stack  is  located  to  the  right  of  the  dotted  line. 
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Resonance  frequency  vs  DeltaT  for  the  annular  prime  mover 
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Figure  4.5  Comparison  of  calculated  and  measured  resonance  frequencies  vs.  AT 
for  the  annular  prime  mover. 
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Figure  4.6    Comparison  of  calculated  and  measured  \IQ  vs.  AT  for  the  annular 
prime  mover. 
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Low  mode  with  dT  =  0  Kelvin  ;  freq  =  423.1+5.735i 
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Figure  4.7(a)  Mode  shape  for  the  low  mode  of  the  annular  prime  mover  when  the 
driver  is  located  90°  from  the  stack  and  AT  =  0  K. 
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Low  mode  with  dT  =  79.63  Kelvin  ;  freq  =  427.4+5.722i 
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Figure  4.7(b)    Mode  shape  for  the  low  mode  of  the  annular  prime  mover  when 
the  driver  is  located  90°  from  the  stack  and  AT  =  80  K. 
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Low  mode  with  dT  =  1 99.6  Kelvin  ;  freq  =  438.3+5.442i 
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Figure  4.7(c)  Mode  shape  for  the  low  mode  of  the  annular  prime  mover  when  the 
driver  is  located  90°  from  the  stack  and  AT  =  200  K. 
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High  mode  with  dT  =  0  Kelvin  ;  freq  =  438.3+2.099i 
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Figure  4.8(a)   Mode  shape  for  the  high  mode  of  the  annular  prime  mover  when 
the  driver  is  located  1 80°  from  the  stack  and  AT  =  0  K. 
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High  mode  with  dT  =  75.73  Kelvin  ;  freq  =  442.7+2.245i 
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Figure  4.8(b)   Mode  shape  for  the  high  mode  of  the  annular  prime  mover  when 
the  driver  is  located  180°  from  the  stack  and  AT  =  76  K. 
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High  mode  with  dT  =  192.1  Kelvin  ;  freq  =  453.5+2.962i 
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Figure  4.8(c)   Mode  shape  for  the  high  mode  of  the  annular  prime  mover  when 
the  driver  is  located  180°  from  the  stack  and  AT  =  192  K. 
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B.     CONSTRICTED  ANNULAR  PRIME  MOVER 

In  this  section,  results  for  the  annular  prime  mover  with  a  constriction  located  90° 
from  the  center  of  the  stack/heat  exchanger  assembly  are  presented.  The  constriction  is  45° 
long  and  three  different  porosities  0.7,  0.3,  and  0.1  are  used.  First,  the  comparisons  of  the 
measured  and  calculated  resonance  frequencies  of  the  constricted  prime  mover  are 
presented.  Next,  the  comparisons  of  the  measured  and  calculated  Q  of  the  constricted 
prime  mover  are  discussed.  The  end  corrections  are  included  in  the  calculations  using  the 
two  methods  discussed  in  Chapter  II.  This  section  concludes  with  the  comparisons  of  the 
measured  and  computed  mode  shapes  for  the  constricted  prime  mover. 

As  mentioned  earlier,  the  measured  resonance  frequencies  and  Q's  are  determined 
from  pole-zero  analysis  of  the  frequency  response.  In  the  measurement,  the  driver  is 
located  45°  from  the  center  of  the  stack  assembly  and  the  constriction.  The  theoretical 
resonance  frequencies  and  Q's  were  calculated  for  four  different  cases,  represented  by  A, 
B,  C,  and  D.  All  the  results  are  listed  in  Tables  4.1  and  4.2.  In  these  tables,  Column  A 
represents  the  calculated  results  without  any  end  corrections.  Columns  B  and  C  represent 
the  results  with  end  corrections  for  the  constriction  using  the  conformal  transformation 
method  and  higher  order  mode  method,  respectively.  It  is  worth  pointing  out  that  there  are 
no  appreciable  differences  between  of  the  results  for  \IQ  calculated  from  the  conformal 
mapping  transformation  and  the  higher  order  mode  method  when  the  end  corrections  in  the 
constriction  are  included.  Column  D  gives  the  results  for  the  case  of  applying  end 
corrections  to  both  the  constriction  and  the  stack  using  the  higher  order  mode  method.  The 
end  corrections  in  the  prime  mover  stack  are  approximated  by  applying  the  end  corrections 
to  one  slit  and  summing  the  impedance  of  the  slits  in  parallel.  The  temperature  differences 
were  determined  by  Eqs.  (4.1)  and  (4.2).    Table  4.1  lists  the  calculated  and  measured 
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resonance  frequencies  for  the  low  and  high  modes  for  the  various  constrictions.  Table  4.2 
lists  the  calculated  and  measured  Q  for  the  low  and  high  modes  for  the  various 
constrictions.  In  these  tables,  Rs  represents  the  area  ratio  of  the  constriction,  i.e.,  the 
porosity. 

It  is  seen  that  the  inclusion  of  end  corrections  has  significant  effects  on  the  predictions 
for  the  low  mode.  For  example,  in  the  case  of  Rs  =0.1,  the  averaged  differences 
between  the  calculated  and  measured  resonance  frequencies  are  1.94  %  (no  end 
corrections),  0.67%  (end  corrections  using  conformal  transformation  method),  and  0.39% 
(end  corrections  using  higher  order  mode  method).  On  the  contrary,  including  end 
corrections  for  the  constriction  has  little  effect  on  the  predictions  for  the  high  mode.  This 
point  is  evident  from  Figs.  4.9  and  4.10,  where  it  is  seen  that  inclusion  of  end  corrections 
in  the  constriction  is  important  for  the  low  mode  only.  However,  inclusions  of  the  end 
corrections  in  both  the  constriction  and  stack  has  more  effect  on  the  high  mode  than  on  the 
low  mode. 

The  point  to  be  taken  from  these  figures  is  that  the  program  does  a  good  job  of 
predicting  the  resonance  frequencies,  even  in  the  absence  of  including  end  corrections.  We 
have  confidence  in  the  end  corrections  applied  to  the  constriction.  This  confidence  is  the 
result  of  preliminary  work  by  Choe  (Choe,  1997).  Application  of  end  corrections  to  the 
stack  is  only  intended  to  be  approximate,  to  indicate  that  they  should  have  an  effect.  The 
contribution  of  this  work  is  to  point  out  the  necessity  of  including  end  corrections.  A 
detailed  analysis  of  stack  end  corrections  is  an  area  for  future  work. 

Figure  4.11  shows  the  comparisons  of  the  calculated  and  measured  1/(2  of  the  low 
mode  for  the  unconstricted  and  constricted  prime  mover  with  various  porosities.  In  this  set 
of  experiments,  when  the  constriction  porosity  is  0. 1 ,  the  prime  mover  reaches  onset  at  a 
temperature  difference  of  260  K.     The  horizontal  error  bar  in  Fig.  4.11   indicates  the 
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uncertainty  in  determining  the  onset  temperature  from  the  measure  temperature  distribution 
as  discussed  previously.  In  all  cases,  the  calculated  results  under  estimate  l/Q.  Therefore, 
the  calculated  results  show  a  significantly  lower  onset  value  of  AT  in  the  case  of  Rs  =  0.1. 
The  reason  for  this  discrepancy  is  not  fully  understood.  Nonetheless,  both  the  measured 
and  calculated  l/Q  indicate  the  same  trends.  First,  the  relative  values  for  l/Q  are  consistent. 
Also,  both  the  predicted  and  measured  results  for  the  constricted  prime  mover  intersect 
within  a  narrow  range  of  AT.  The  predicted  values  intersect  at  AT  =  100K  and  the 
measured  data  intersect  at  AT  =  175  K.  The  l/Q  curve  for  Rs  =0.1  begins  with  the 
highest  value  of  the  three  constricted  cases  and  ends  with  the  lowest  value  after  the 
intersection  at  AT  =100  K.  In  contrast,  l/Q  for  Rs  =  0.7  begins  with  the  lowest  value 
and  ends  with  the  largest  value  at  higher  AT.  l/Q  for  Rs  =  0.3  remains  in  between  l/Q 
for  the  other  two  porosities. 

The  conclusion  to  be  drawn  from  Fig.  4.11  is  that  it  is  predicted  that  a  constricted 
annular  prime  mover  will  reach  onset,  and  it  does.  Although,  the  agreement  between  the 
measured  and  calculated  results  are  only  fair,  there  is  very  good  qualitative  agreement. 
This  indicates  that  although  there  are  still  some  details  to  be  investigated,  the  general 
behavior  is  understood. 

Figure  4.12  shows  the  comparisons  between  the  calculated  and  measured  l/Q  for  the 
high  mode  for  the  unconstricted  and  constricted  prime  mover.  It  is  seen  that  the 
temperature  difference  has  less  effect  on  l/Q  for  the  high  mode  than  for  the  low  mode. 
Over  the  temperature  span  from  AT  =0  K  to  AT  =  200  K,  the  change  in  l/Q  is  less  than 
30  %  in  all  cases.  The  same  comments  about  the  general  qualitative  agreement  apply  to  the 
high  mode  results  as  they  do  to  the  low  mode  results. 

The  question  to  be  addressed  now  is  why  does  a  constricted  annular  prime  mover 
reach  onset.  The  answer  lies  in  the  effects  that  a  constriction  has  on  mode  shapes. 


93 


Figures  4.13(a)  and  (b)  show  the  comparison  of  the  measured  and  computed  mode 
shapes  of  the  low  mode  for  Rs  =  0.1.  In  this  case  the  constriction  is  much  longer  than  the 
stack  assembly  (LconstrictlJLslack=  3.33)  and  the  porosity  of  the  constriction  is  much  smaller 
than  that  of  the  stack  assembly.  As  a  result,  the  mode  shape  is  dominated  by  the 
constriction.  There  is  a  pressure  node  near  the  center  of  the  constriction  and  a  pressure 
antinode  between  the  constriction  and  the  stack  assembly  nearest  (but  not  at)  the  hot  end  of 
the  stack.  End  corrections  for  the  constriction  are  included  using  the  higher  order  mode 
method.  Figure  4.13(a)  is  for  AT  =  0  K  and  Fig.  4.13(b)  is  for  AT  =288  K.  The 
frequencies  are  the  calculated  values.  It  is  observed  that  there  is  a  good  agreement  between 
the  measured  and  predicted  mode  shape  at  AT  =0  K.  This  mode  shape  is  favorable  to 
thermoacoustic  growth.  The  agreement  in  mode  shape  worsens  at  AT  =288  K.  It  should 
be  noticed  that  the  imaginary  part  of  the  calculated  resonance  frequency  is  negative,  which 
suggests  the  onset  has  been  reached. 

Figures  4.14(a)  and  (b)  show  the  comparison  of  the  measured  and  computed  mode 
shapes  of  the  high  mode  for  Rs  =0.1.  The  agreement  between  the  measured  and 
calculated  mode  shape  is  good.  As  one  would  expect,  there  is  a  pressure  antinode  near  the 
center  of  the  constriction. 

One  important  common  feature  of  Figs.  4.13  and  4.14  is  that  the  mode  shapes  are 
essentially  determined  by  the  constriction  and  are  independent  of  the  temperature 
difference.  This  is  a  key  to  making  the  annular  prime  mover  reach  onset. 

The  main  conclusion  from  this  section  is  that  to  build  a  functional  annular  prime 
mover,  the  relative  locations  of  the  constriction  and  the  stack,  and  their  relative  porosities 
are  the  crucial  factors.  This  is  equivalent  to  incorporating  some  type  of  dominating 
boundary  conditions  into  the  prime  mover  to  alter  the  mode  shapes. 
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418.2 

77.1 

421.9 

421 

420.5 

420.4 

418.8 

112.2 

423.5 

422.5 

422.1 

421.9 

419.9 

179.2 

426.9 

425.9 

425.4 

425.2 

422 

231 

433.9 

429.6 

429.2 

428.8 

424.5 

Rs=0.7   High   mode 

A 

B 

C 

D 

DeltaT  (K) 

f,  cal  (Hz) 

f,  cal  (Hz) 

f,  cal  (Hz) 

f,  cal  (Hz) 

f,  mea  (Hz) 

0 

438.9 

438.8 

438.7 

434.2 

436.2 

46.8 

442.5 

442.4 

442.3 

438.1 

440.2 

77.1 

445.2 

445.1 

445.0 

440.9 

442.2 

112.2 

448.2 

448.1 

448.0 

444 

444.4 

179.2 

454.2 

454.1 

454.0 

450.2 

448.5 

231 

460.2 

460 

460.0 

456.3 

452 

Table  4.1  Comparisons  of  the  calculated  and  measured  resonance  frequencies  for  the 
constricted  annular  prime  mover.  Case  A:  without  end  corrections;  Cases  B,  C,  D: 
End  corrections  with  different  methods. 
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Rs=0.1    Low    mode 

A 

B 

C 

D 

DeltaT  (K) 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qlow,  mea 

0 

50.4 

45.5 

45.4 

38.2 

43.3 

41.1 

68.8 

60.8 

60.9 

48.9 

51.5 

70.8 

92.3 

80.2 

80.2 

61.3 

60.5 

103.9 

148.3 

122.6 

123.4 

84.8 

73.8 

175.4 

-519.6 

-1089 

-1029 

421 

127 

228.2 

-124.3 

-123.3 

-131.9 

-225 

239.5 

Rs=0.1    High   mode 

A 

B 

C 

D 

DeltaT  (K) 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qhigh,  mea 

0 

36.4 

36.4 

36.4 

24.3 

32.6 

41.1 

37.9 

37.9 

37.9 

24.9 

33.7 

70.8 

39.2 

39.1 

39.1 

25.4 

34.2 

103.9 

40.8 

40.8 

40.8 

26 

35.9 

175.4 

45.7 

45.6 

45.7 

27.8 

37.9 

228.2 

51.3 

51.2 

51.2 

29.6 

41.3 

Rs=0.3   Low   mode 

A 

B 

C 

D 

DeltaT  (K) 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qlow,  mea 

0 

74.9 

69.6 

69.1 

62.7 

65.1 

42.1 

90.2 

84.8 

83.7 

75.6 

73.8 

71.7 

105.5 

98.6 

98.1 

87.7 

81.1 

105.6 

129.2 

122 

123 

107 

90.9 

175.4 

235.2 

234 

248 

196 

121.7 

226.4 

514.4 

590 

690 

435 

157 

Rs=0.3    High   mode 

A 

B 

C 

D 

DeltaT  (K) 

Qhigh,  cal 

Qhigh,  cal 

Qhigh,  cal 

Qhigh,  cal 

Qhigh,  mea 

0 

35.7 

35.7 

35.7 

24 

32.8 

42.1 

36.9 

36.9 

36.9 

24.5 

33.7 

71.7 

37.9 

37.8 

37.8 

24.9 

34.3 

105.6 

39.2 

39.1 

39.1 

25.3 

35.1 

175.4 

42.6 

42.4 

42.4 

26.6 

37.4 

226.4 

45.8 

45.8 

45.8 

27.8 

38.8 

Rs=0.7  Low   mode 

A 

B 

C 

D 

DeltaT  (K) 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qlow,  cal 

Qlow,  mea 

0 

98.8 

97.3 

97.2 

94.2 

85.4 

46.8 

108.2 

106.8 

106.6 

103.4 

95 

77.1 

118.4 

116.8 

116.7 

114.4 

98.2 

112.2 

136.7 

134.5 

134.5 

133.6 

103 

179.2 

216.3 

207 

205.8 

210 

114 

231 

312.5 

407.1 

398.2 

444.3 

119.1 

Rs=0.7   High    mode 

A 

B 

C 

D 

DeltaT  (K) 

Qhigh,  cal 

Qhigh,  cal 

Qhigh,  cal 

Qhigh,  cal 

Qhigh,  mea 

0 

36.6 

36.6 

36.6 

24.7 

35.2 

46.8 

36.5 

36.5 

36.5 

24.7 

36.4 

77.1 

36.1 

36.1 

36.2 

24.3 

36.3 

112.2 

35.7 

35.7 

35.8 

24 

35.2 

179.2 

34.5 

34.8 

34.9 

23.4 

36.4 

231 

33.7 

34 

34.1 

22.8 

37.2 

Table  4.2  Comparisons  of  the  calculated  and  measured  Q  for  the  constricted 
annular  prime  mover.  Case  A:  without  end  corrections;  Cases  B,  C,  D:  End 
corrections  with  different  methods. 
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315 


Low  mode,  Resonance  frequency  vs.  DeltaT ,  Rs  =  0.1 
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Calculated,  no  end  corrections 

Calculated,  end  corrections  at  constriction(A) 

Calculated,  end  corrections  at  constriction(B) 

-  -  Calculated,  end  corrections  at  constriction  &  stack  (B) 
o     Measured  data 
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Figure  4.9  Comparisons  of  the  measured  and  calculated  resonance  frequencies  of 
the  low  mode  for  the  constricted  prime  mover  with  a  porosity  of  0. 1 .  Method  A 
is  the  conformal  mapping  transformation  and  Method  B  is  the  higher  order  mode. 
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High  mode,  Resonance  frequency  vs.  DeltaT  ,  Rs  =  0.1 
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Calculated,  no  end  corrections 

Calculated,  end  corrections  at  constriction(A) 

Calculated,  end  corrections  at  constriction(B) 

-  -  Calculated,  end  corrections  at  constriction  &  stack  (B) 
o     Measured  data 


0 


50 


100 

Delta  T  (Kelvin) 


150 


200 


Figure  4.10  Comparisons  of  the  measured  and  calculated  resonance  frequencies 
of  the  high  mode  for  the  constricted  prime  mover  with  a  porosity  of  0. 1 .  Method 
A  is  the  conformal  mapping  transformation  and  Method  B  is  the  higher  order 
mode. 
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1/Q  for  the  low  mode  of  the  annular  prime  mover 
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Figure  4. 1 1  Comparisons  of  the  measured  and  predicted  \IQ  of  the  low  mode  for 
the  unconstricted  and  constricted  prime  mover.      Symbol   o    represents   the 

measured  \IQ  of  unconstricted  prime  mover.  Symbols  *,  +,  and  x  represent  the 

measured  \IQ  for  the  constricted  prime  mover  with  porosities  of  0.1,  0.3,  and 
0.7,  respectively.  The  lines  represent  the  predicted  values.  No  end  corrections 
are  included  in  the  predicted  values  for  the  unconstricted  prime  mover. 
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1/Q  for  the  high  mode  of  the  annular  prime  mover 
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Figure  4.12  Comparisons  of  the  measured  and  predicted  l/Q  of  the  high  mode 
for  the  unconstricted  and  constricted  prime  mover.     Symbol  o  represents  the 

measured  l/Q  of  unconstricted  prime  mover.  Symbols  *,  +,  and  x  represent  the 

measured  l/Q  for  the  constricted  prime  mover  with  porosities  of  0.1,  0.3,  and 
0.7,  respectively.  The  lines  represent  the  predicted  values.  No  end  corrections 
are  included  in  the  predicted  values  for  the  unconstricted  prime  mover. 
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Low  mode  with  end  effect  2  ;  dT  =  0  Kelvin  ;  freq  =  296.2+3.259i ;  Rs  =  0.1 


0.4 
Position,  m 


Figure  4.13(a)    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs  =  0.1)  when  the  driver  is  located  45°  from  the  stack  and  AT  =  0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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Low  mode  with  end  corrections  ;  dT  =  228.2  Kelvin  ;  freq  =  307.4-1.165i ;  Rs  =  0.1 
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Figure  4.13(b)    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs  =  0.1)    when  the  driver  is  located  45°  from  the  stack  and  AT  =  228  K. 

The  calculated  results  are  based  on  the  higher  order  mode  method.  The  frequency 
is  the  calculated  value. 
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High  mode  with  end  effect ;  dT  =  0  Kelvin  ;  freq  =  474+6.51 4i ;  Rs  =  0.1 
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0.3  0.4 

Position,  m 

Figure  4.14(a)    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs  =  0.1)  when  the  driver  is  located  45°  from  the  stack  and  AT  =  0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 


103 


High  mode  with  end  effect  2  ;  dT  =  228.2  Kelvin  ;  freq  =  493.1+4.817i ;  Rs  =  0.1 


0.3  0.4 
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Figure  4. 14(b)    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs  =  0.1)    when  the  driver  is  located  45°  from  the  stack  and  AT  =  228  K. 

The  calculated  results  are  based  on  the  higher  order  mode  method.  The  frequency 
is  the  calculated  value. 
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C.    TWO  STACK  ANNULAR  PRIME  MOVER 

The  knowledge  gained  from  the  last  two  sections  leads  to  the  design  of  a  symmetric 
two  stack  annular  prime  mover  in  which  the  constriction  in  the  last  section  is  replaced  by  a 
stack  assembly.  This  proposed  design  of  the  annular  prime  mover,  as  shown  in  Fig.  4. 15, 
has  two  identical  stack/heat  exchanger  assemblies.  The  two  stack  assemblies  are  placed 
such  that  the  hot  heat  exchangers  face  each  other. 

The  MATLAB  program  was  modified  to  predict  the  performance  of  this  two  stack 
annular  prime  mover,  by  adding  one  more  guess  of  the  time-averaged  energy  flux  Hi 
along  the  second  stack  to  match  the  temperature  of  the  hot  heat  exchanger.  Two  different 
cases  are  studied  here.  In  case  1  the  duct  between  the  two  heat  exchangers  is  held  at 
ambient  temperature  Ta  =  293  K  while  in  case  2  it  is  held  at  the  desired  temperature  of  the 
hot  heat  exchanger  Th.  The  reason  for  studying  case  1  is  that  the  current  apparatus  can  be 
easily  modified  to  conduct  the  experiment.  However,  because  of  the  temperature 
discontinuity  near  the  junction  of  the  duct  and  hot  heat  exchanger,  it  is  difficult  to  model 
this  case  theoretically  and  computationally.  It  is  known  that  the  temperature  discontinuity 
in  the  duct  will  also  create  a  small  discontinuity  in  the  work  flow  (Rott,  1969).  On  the 
contrary,  although  the  theoretical  analysis  of  case  2  is  simpler  than  case  1,  investigating  it 
requires  a  new  apparatus. 

Figure  4.16  shows  the  calculated  resonance  frequencies  of  the  two-stack  prime  mover 
for  the  two  cases.  In  case  1,  it  is  seen  that  the  resonance  frequencies  for  the  high  mode 
increase  approximately  linearly  with  temperature  difference.  However,  AT  has  little  effect 
on  the  resonance  frequencies  of  the  low  mode  in  case  1 .  The  resonance  frequencies  for 
both  the  low  and  high  mode  in  case  2  increase  approximately  linearly  with  temperature 
difference.    Because  the  duct  between  the  two  hot  heat  exchangers  is  held  at  7h,  the 
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increase  in  resonance  frequency  is  bigger  than  that  in  case  1 .  It  should  be  pointed  out  that 
in  case  1  at  very  low  temperature  difference  (AT  <  50  K),  the  resonance  frequencies  of 
the  low  mode  and  high  mode  are  very  nearly  equal.  The  MATLAB  program  always 
converges  to  only  one  mode.  For  example,  at  AT  =  0  K,  the  program  always  converges 
to  the  high  mode,  while  at  AT  =  25  K  it  always  converges  to  the  low  mode.  As  AT 
increases,  the  modes  are  well  separated  and  no  such  condition  was  observed.  This 
condition  is  not  present  in  case  2.  It  was  always  possible  to  distinguish  between  the  low 
mode  and  high  mode.  It  should  also  be  pointed  out  that  for  case  1  at  AT  =  0  K,  the 
resonance  frequency  for  the  high  mode  is  424.6  Hz  and  is  424.0  Hz  for  the  low  mode. 
The  results  indicate  the  possibility  of  mode  level  repulsion  (Pippard,  1985  and  1989; 
Laraza  and  Denardo  1997). 

Figure  4.17  shows  the  calculated  \IQ  of  the  low  mode  and  high  mode  for  the  two- 
stack  prime  mover  for  the  two  cases.  It  is  evident  that  the  high  mode  reaches  onset  at 
AT  =  190  K  for  case  1.  The  onset  temperature  difference  for  the  high  mode  in  case  2  is 
approximately  350  K.  An  interesting  feature  different  from  the  \IQ  for  the  constricted 
prime  mover  is  the  linearly  increase  and  decrease  in  \IQ  for  the  low  mode  and  high  mode, 
respectively.  Similar  dependence  has  been  observed  for  a  rigid-rigid  prime  mover  (Lin, 
1989;  Atchley,  et  al.,  1992). 

Frequency  splitting  in  the  two-stack  annular  prime  mover  offers  the  possibility  that  the 
harmonics  of  the  fundamental  frequency  and  the  overtones  of  the  resonator  will  not 
coincide,  thus  detuning  the  system.  This  means  that  if  high  amplitudes  are  achieved  above 
onset,  waveform  distortion  and  shocking  should  not  be  a  problem  (Coppens,  1971; 
Atchley  and  Hofler,  1990;  Atchley  et.  al.  1993). 

Figures  4.18  to  4.21  show  the  calculated  mode  shapes  for  the  low  mode  and  high 
mode  for  the  two  cases.  As  one  would  expect  from  symmetry,  there  is  a  pressure  antinode 
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between  the  two  stacks  for  the  high  mode  and  a  pressure  node  between  the  two  stack 
assembly  for  the  low  mode.  It  is  evident  that  only  the  high  modes  in  both  cases  support 
thermoacoustic  growth,  as  is  confirmed  by  Fig.  4.17. 


Hot  Heat 
Exchanger  (Th) 


Prime  Mover 
Stack 


Ambient  Heat 
Exchanger  (Ta) 


Prime  Mover 
Stack 


Ambient  Heat 
Exchanger  (Ta) 


Figure  4. 15  The  two-stack  annular  thermoacoustic  prime  mover. 
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Calculated  resonance  frequency  vs.  dT  for  the  two  stack  annular  prime  mover 
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Figure  4.16  The  calculated  resonance  frequencies  for  the  two  stack  annular  prime 
mover.  Case  1:  The  duct  between  the  two  hot  heat  exchangers  is  held  at  Ta. 
Case  2:  The  duct  between  the  two  hot  heat  exchangers  is  held  at  Th. 
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Calculated  1/Q  vs.  dT  for  the  two  stack  annular  prime  mover 
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Figure  4.17  The  calculated  \IQ  for  the  low  mode  and  high  mode  for  the  two 
stack  annular  prime  mover.  Case  1:  The  duct  between  the  two  hot  heat 
exchangers  is  held  at  Ta.  Case  2:  The  duct  between  the  two  hot  heat  exchangers 
is  held  at  Th. 
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Mode  shape  for  of  the  two  stack  prime  mover,  dT  =  200K;  freq  =  426.1+14.74i 
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Figure  4.18  The  calculated  mode  shape  of  the  low  mode  of  the  two  stack  annular 
prime  mover  at  AT  =  200  K.  Case  1:  The  duct  between  the  two  hot  heat 
exchangers  is  held  at  Ta  =  293  K. 
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Mode  shape  for  of  the  two  stack  prime  mover,  dT  =  199.2K;  freq  =  433.5-0.61 4i 
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Figure  4.19  The  calculated  mode  shape  of  the  high  mode  of  the  two  stack  annular 
prime  mover  at  AT  =  200  K.  Case  1:  The  duct  between  the  two  hot  heat 
exchangers  is  held  at  Ta  =  293  K. 
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Mode  shape  for  of  the  two  stack  prime  mover,  dT  =  200K;  freq  =  535.9+1 8.2i 
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Figure  4.20  The  calculated  mode  shape  of  the  low  mode  of  the  two  stack  annular 
prime  mover  at  AT  =  200  K.  Case  2:  The  duct  between  the  two  hot  heat 
exchangers  is  held  at  Th  =  493  K. 
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Mode  shape  for  of  the  two  stack  prime  mover,  dT  =  199.9K;  freq  =  540.7+3.941  i 
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Figure  4.2 1  The  calculated  mode  shape  of  the  high  mode  of  the  two  stack  annular 
prime  mover  at  AT  =  200  K.  Case  2:  The  duct  between  the  two  hot  heat 
exchangers  is  held  at  Th  =  493  K. 
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V.     SUMMARY  AND  CONCLUSION 

The  objective  of  this  dissertation  is  to  investigate  a  thermoacoustic  prime  mover  having 
periodic  boundary  conditions  both  experimentally  and  theoretically.  The  experimental 
approach  has  been  to  design,  build,  and  test  an  annular  thermoacoustic  prime  mover.  The 
resonance  frequencies,  quality  factors,  and  the  mode  shapes  of  the  prime  mover  are 
measured  as  functions  of  the  temperature  difference  across  the  prime  mover  stack.  The 
experimental  resonance  frequencies  and  Q's  are  determined  from  pole-zero  analysis  of  the 
measured  frequency  response. 

An  important  aspect  of  this  dissertation  is  the  computational  analysis  of  a 
thermoacoustic  device  with  periodic  boundary  conditions.  Two  approaches  were  taken. 
First  an  existing  program,  DeltaE,  was  applied  to  this  novel  geometry.  It  was  found  that 
DeltaE  can  be  applied  to  the  problem,  but  has  some  disadvantages.  It  is  difficult  to  extract 
Q's  below  onset  and  awkward  to  extract  mode  shapes.  It  is  also  difficult  to  include 
temperature  gradients  in  ducts. 

Because  of  these  problems,  a  new  program  was  developed.  It  uses  the  same  basic 
approach  as  DeltaE,  but  overcomes  the  disadvantages.  The  complex  eigenfrequency  is 
used  to  extract  resonance  frequencies  and  Q  's. 

The  results  of  this  dissertation  indicate  that  under  all  but  perhaps  extreme  conditions  an 
annular  prime  mover  will  not  reach  onset.  The  reason  is  well  understood.  The  eigenmodes 
of  a  single-stack  annular  prime  mover  do  not  support  thermoacoustic  growth.  This  finding 
lead  to  the  investigation  of  a  constricted  annular  prime  mover.  The  idea  is  that  the 
constriction  could  impose  dominating  boundary  conditions  on  the  prime  mover  thus 
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altering  the  eigenmodes.  It  was  shown  that  decreasing  the  porosity  of  the  constriction 
forced  the  prime  mover  into  onset.  The  porosity  and  location  of  the  constriction  is  used  to 
adjust  the  eigenmodes  of  the  system  so  that  thermoacoustic  growth  is  supported. 

Another  way  to  produce  eigenmodes  that  support  thermoacoustic  growth  is  to 
introduce  a  symmetry  into  the  system  which  shifts  the  nodes  and  antinodes  of  the  modes 
off  the  stacks,  such  that  for  one  of  the  modes  a  pressure  antinode  is  near  the  hot  end  of  the 
stacks.  It  is  predicted  that  a  two-stack  annular  prime  mover  will  reach  onset.  This  design 
is  an  area  for  future  research. 

Another  area  for  future  study  is  end  correction  for  stacks.  It  was  shown  in  this 
dissertation  that  inclusion  of  end  corrections  does  affect  the  predicted  resonance  frequencies 
and  <2's.  Determining  the  proper  end  corrections  for  realistic  stack  geometries  is  likely  to 
be  a  nontrivial  project. 
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APPENDIX  A.l  THE  DELTAE  INPUT  FILE  FOR  THE  LOW  FREQUENCY 
MODE  OF  THE  ANNULAR  THERMOACOUSTIC  PRIME  MOVER 


TITLE  Annular  resonator   (file  :  to2dTflowC.in) 
!Created@23: 18:35  02-MAY-97  with  DeltaE  Vers.  3 
! —    0 

BEGIN       Initial 
1.0130E+05       a  Mean  P     Pa 

b  Freq.       Hz 

c  T-beg      K 

d  lpl@0      Pa 

e  Ph(p)0    deg 

f  IUI@0     mA3/s 

gPh(U)0    deg 
Gas  type 

Solid  type 


Obi  for  the  Power  Macintosh 


427.0 

293.0 

16.63 

-43.90 

6.7062E-04 

38.50 

air 

ideal 

t „ 

SOFTEND 

0.0000 

0.0000 


16.63 

A  lpl@0  G(  Od) 

I 

-43.90 

B  Ph(p)0  G(  Oe) 

P 

6.7062E-04 

CIUI@0  G(0f) 

P 

G 

38.50 

D  Ph(U)0  G(  Og) 

P 

G 

-1.8309E-03 

E  Heatln  G(27e) 

P 

G 

G 

1 


1 


a  Re(Z) 
b  Im(Z) 


(t) 
(t) 


Gas  type 
Solid  type 


sameas  0 

ideal 

i 

ISODUCT     Duct 
2.6300E-03       a  Area 
0.2050 
3.1788E-02 


mA2 
b  Perim       m 
c  Length       m 


16.63 
-43.90 
6.7062E-04 
38.50 

7.3673E-04 
7.3673E-04 
2.0847E-02 
-0.1564 
293.0 


sameas  0 


Gas  type 


42.05 
-48.73 

6.2404E-04 

38.20 

7.0088E-04 


A  Ipl         Pa 
B  Ph(p)       deg 
C  IUI       mA3/s 
D  Ph(U)       deg 
EHdot 
FWork 
G  Re(Z) 
H  Im(Z) 
I    T        K 


W 
W 


A  Ipl         Pa 
B  Ph(p)      deg 
C  IUI      mA3/s 
D  Ph(U)       deg 
E  Hdot       W 
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ideal 

i __ 

Solid  type 

...       7,    ... 

7.0088E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

64.95 

Alpl 

Pa 

0.2050 

b  Perim 

m 

-49.96 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

5.3889E-04 

CIUI      i 

mA3/s 

37.86 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

6.6858E-04 

EHdot 

W 

ideal 

i 

Solid  type 

___      A    ... 

6.6858E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

83.85 

Alpl 

Pa 

0.2050 

b  Perim 

m 

-50.55 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

4.2045E-04 

CIUI      i 

mA3/s 

37.36 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

6.4101E-04 

EHdot 

W 

ideal 

i 

Solid  type 

<N      ... 

6.4101E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

97.57 

Alpl 

Pa 

0.2050 

b  Perim 

m 

-50.93 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

2.7606E-04 

CIUI 

mA3/s 

36.44 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

6.1820E-04 

EHdot 

W 

ideal 
t 

Solid  type 

...     6  - 

6.1820E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

105.3 

Alpl 

Pa 

0.2050 

b  Perim 

m 

-51.21 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

1.1484E-04 

CIUI 

mA3/s 

33.10 

D  Ph(U) 

deg 

sameas  0  Gas  type 

5.9903E-04 

EHdot 

W 

ideal       Solid  type 

! 

. —      7    .. 

5.9903E-04 

FWork 

W 

ISODUCT     Duct 
2.6300E-03       a  Area 


mA2 


106.4 


A  Ipl        Pa       P 
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0.2050 
3.1788E-02 


b  Perim       m 
c  Length       m 


sameas  0  Gas  type 

ideal  Solid  type 

! 8 

ISODUCT  Duct 

2.6300E-03  a  Area       mA2 

0.2050  b  Perim       m 

3.1788E-02  c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

i 9 

ISODUCT     Duct  7 

2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02      c  Length      m 


-51.45 

B  Ph(p) 

deg 

5.5371E-05 

CIUI      i 

mA3/s 

-130.1 

D  Ph(U) 

deg 

5.8148E-04 

EHdot 

W 

5.8148E-04 

FWork 

W 

101.0 

Alpl 

Pa 

-51.67 

B  Ph(p) 

deg 

2.1977E-04 

CIUI 

mA3/s 

-138.8 

D  Ph(U) 

deg 

5.6316E-04 

EHdot 

W 

5.6316E-04 

FWork 

W 

sameas  0 

Gas  type 

5.4183E-04 

EHdot 

W 

ideal 

i 

Solid  type 
in  — 

5.4183E-04 

FWork 

W 

ISODUCT 

Duct               8 

2.6300E-03 

a  Area       mA2 

72.22 

Alpl 

Pa 

0.2050 

b  Perim       m 

-52.23 

B  Ph(p) 

deg 

3.1788E-02 

c  Length      m 

4.9997E-04 

CIUI 

mA3/s 

-140.6 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

5.1603E-04 

EHdot 

W 

ideal 

i _ 

Solid  type 
1 1   __ 

5.1603E-04 

FWork 

w 

ISODUCT 

Duct 

2.6300E-03 

a  Area       mA2 

50.59 

Alpl 

Pa 

0.2050 

b  Perim       m 

-52.75 

B  Ph(p) 

deg 

3.1788E-02 

c  Length      m 

5.9773E-04 

CIUI 

mA3/s 

-140.9 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

4.8535E-04 

EHdot 

W 
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ideal 

i 

Solid  type 

—   12  — 

4.8535E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

25.84 

Alpl 

Pa 

0.2050 

b  Perim 

m 

-54.17 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

6.5855E-04 

CIUI      i 

mA3/s 

-141.1 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

4.5059E-04 

EHdot 

W 

ideal 

i 

Solid  type 

.„  ii  ... 

4.5059E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

1.340 

Alpl 

Pa 

0.2050 

b  Perim 

m 

-165.9 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

6.7865E-04 

CIUI 

mA3/s 

-141.3 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

4.1353E-04 

EHdot 

W 

ideal 

! 

Solid  type 

...  14  ... 

4.1353E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

26.91 

Alpl 

Pa 

0.2050 

b  Perim 

m 

131.0 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

6.5679E-04 

CIUI 

mA3/s 

-141.5 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.7652E-04 

EHdot 

W 

ideal 

t 

Solid  type 

....  15  .. 

3.7652E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

51.56 

Alpl 

Pa 

0.2050 

b  Perim 

m 

129.7 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

5.9432E-04 

CIUI 

mA3/s 

-141.6 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.4190E-04 

EHdot 

W 

ideal 

i 

Solid  type 

—  16  - 

3.4190E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

73.03 

Alpl 

Pa 
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0.2050 

b  Perim 

m 

129.2 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

4.951 1E-04 

CIUI       i 

T1A3/S 

-141.8 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.1142E-04 

EHdot 

W 

ideal 

i 

Solid  type 

....  n  ... 

3.1142E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

89.99 

Alpl 

Pa 

0.2050 

b  Perim 

m 

128.9 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

3.6528E-04 

CIUI      i 

mA3/s 

-142.1 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

2.8583E-04 

EHdot 

W 

ideal 

i 

Solid  type 

....  is  ... 

2.8583E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

101.4 

Alpl 

Pa 

0.2050 

b  Perim 

m 

128.8 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

2.1289E-04 

CIUI      i 

mA3/s 

-142.6 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

2.6468E-04 

EHdot 

W 

ideal 

i 

Solid  type 

....  19  ... 

2.6468E-04 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

106.5 

Alpl 

Pa 

0.2050 

b  Perim 

m 

128.7 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

4.7479E-05 

CIUI 

mA3/s 

-146.9 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

2.4644E-04 

EHdot 

W 

ideal 

i 

Solid  type 

— -  20  — 

2.4644E-04 

FWork 

W 

VDUCER 

Driver 

1.0000E-09 

a  Re(Ze) 

ohms 

106.5 

Alpl 

Pa       P 

0.0000 

b  Im(Ze) 

128.7 

B  Ph(p) 

deg 

1.0000E+04 

c  Re(Tl)  V-s/mA3 

1.4217E-04 

CIUI 

mA3/s 

0.0000 

d  Im(Tl)  V-s/mA3 

-169.5 

D  Ph(U) 

deg 

-1.0000E+04     eRe(T2) 

Pa/A 

3.5727E-03 

EHdot 

W 
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0.0000 

f  Im(T2) 

Pa/A 

3.5727E-03 

FWork 

w 

1  .OOOOE-09 

g  Re(Zm)  Pa-s/mA3 

3.3263E-03 

G  Workln     W 

1.0000E-09 

h  Im(Zm)  Pa-s/mA3 

1.000 

H  Volts 

V 

1.000 

i  AplVol 

V 

1.0651E-02 

I  Amps 

V 

-51.35 

J  Ph(Ze) 

deg 

sameas  0 

Gas  type 

1.0000E-04 

K  lUxl 

mA3/s 

ideal 

i._ 

Solid  type 

....  21  — - 

-180.0 

L  Ph(-U> 

:    deg 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

108.1 

Alpl 

Pa 

0.2050 

b  Perim 

m 

127.3 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

7.9112E-05 

CIUI      i 

mA3/s 

93.49 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.5543E-03 

EHdot 

W 

ideal 

i 

Solid  type 

.—  22  — ■ 

3.5543E-03 

FWork 

W 

SODUCT     ] 

Duct 

2.6300E-03 

a  Area 

mA2 

103.1 

Alpl 

Pa 

0.2050 

b  Perim 

m 

125.9 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

2.2143E-04 

CIUI      i 

mA3/s 

53.91 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.5352E-03 

EHdot 

W 

ideal 

i 

Solid  type 

....  23  — ■ 

3.5352E-03 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

91.74 

Alpl 

Pa 

0.2050 

b  Perim 

m 

124.2 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

3.7161E-04 

CIUI 

mA3/s 

46.13 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.5132E-03 

EHdot 

W 

ideal 

i 

Solid  type 

....  24  — 

3.5132E-03 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

74.83 

Alpl 

Pa 

0.2050 

b  Perim 

m 

122.0 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

5.0174E-04 

CIUI 

mA3/s 
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42.72 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.4867E-03 

EHdot 

W 

ideal 

i 

Solid  type 
25  — 

3.4867E-03 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area       mA2 

53.48 

Alpl 

Pa 

0.2050 

b  Perim       m 

118.2 

B  Ph(p) 

deg 

3.1788E-02 

c  Length       m 

6.0191E-04 

CIUI      i 

nA3/s 

40.63 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.4553E-03 

EHdot 

W 

ideal 

!_ 

Solid  type 
26  — 

3.4553E-03 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area       mA2 

29.40 

Alpl 

Pa 

0.2050 

b  Perim       m 

108.6 

B  Ph(p) 

deg 

3.1788E-02 

c  Length       m 

6.6543E-04 

CIUI      i 

mA3/s 

39.08 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

3.4195E-03 

EHdot 

W 

ideal 

i 

Solid  type 
27  — 

3.4195E-03 

FWork 

w 

HXFRST 

Ambient  HX 

sameas  3a 

a  Area       mA2 

23.02 

Alpl 

Pa 

0.6290 

b  GasA/A 

103.2 

B  Ph(p) 

deg 

5.0800E-03 

c  Length      m 

6.6930E-04 

CIUI 

mA3/s 

1.7018E-03 

d  yO          m 

38.94 

D  Ph(U) 

deg 

-1.8309E-03 

e  Heatln      W      G 

1.5886E-03 

EHdot 

W 

293.0 

f  Est-T       K      (t) 

3.3455E-03 

FWork 

W 

sameas  0 

Gas  type 

-1.8309E-03 

GHeat 

W 

copper 

Solid  type 

293.0 

H  MetalT      K 

i 

28  — 

STKSLAB 

Prime  Mover  Stack 

2.6300E-03 

a  Area       mA2 

16.26 

Alpl 

Pa 

0.6220 

b  GasA/A 

-43.67 

B  Ph(p) 

deg 

2.4000E-02 

c  Length      m 

6.7080E-04 

CIUI 

mA3/s 

3.0480E-04 

d  yO          m 

38.51 

D  Ph(U) 

deg 

7.5000E-05 

e  Lplate      m 

1.5886E-03 

EHdot 

W 
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7.4220E-04 

FWork 

W 

293.0 

G  T-beg 

K 

sameas  0 

Gas  type 

293.0 

H  T-end 

K 

kapton 

i 

Solid  type 
— 29  

-2.6033E-03 

I  StkWrk 

W 

HXLAST      Cold  HX 

sameas  3a 

a  Area       mA2 

16.63 

A  Ipl         ] 

Pa 

0.7050 

b  GasA/A 

-43.90 

B  Ph(p) 

deg 

3.0480E-04 

c  Length       m 

6.7062E-04 

C  IUI      mA3/s 

1 .2540E-03 

d  yO          m 

38.51 

D  Ph(U) 

deg 

0.0000 

e  Heatln      W      (t) 

7.3668E-04 

EHdot 

W 

293.0 

fEst-T       K    =29H? 

7.3668E-04 

FWork 

W 

sameas  0 

Gas  type 

-8.5195E-04 

GHeat 

W 

copper 

l    

Solid  type 
30  

293.0 

H  MetalT 

K 

SOFTEND 

0.0000 

a  Re(Z)              (t) 

16.63 

A  Ipl 

Pa 

0.0000 

b  Im(Z)              (t) 

-43.90 

B  Ph(p) 

deg 

6.7062E-04 

C  IUI      mA3/s 

38.51 

D  Ph(U) 

deg 

7.3668E-04 

EHdot 

W 

7.3668E-04 

FWork 

W 

2.0846E-02 

G  Re(Z) 

sameas  0 

Gas  type 

-0.1564 

H  Im(Z) 

ideal 

Solid  type 

293.0 

I    T        K 

i 

qi   

DIFFTAR 

Ipl    mismatch 

0.0000 

a  TargDi          =31  A? 

8.0O65E-05 

A  D1-D2 

1A  b  DIAddr 

30A  c  D2Addr 

! _ 

32  — - 

DIFFTAR 

Ph(p)  mismatch 

0.0000 

a  TargDi             =32A?  3.3780E-04 

A  D1-D2 

IB  b  DIAddr 

30B  c  D2Addr 

!._ 

—  33  — - 
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DIFFTAR  IUI    mismatch 

0.0000  aTargDi  =33A?      1.3884E-09    A  D1-D2 

1C  b  DIAddr 

30C  c  D2Addr 

! - — —  34  - - 

DIFFTAR  Ph(U)  mismatch 

0.0000  aTargDi  =34A?     -1.0691E-04  A  D1-D2 

ID  b  DIAddr 

30D  c  D2Addr 

!  The  restart  information  below  was  generated  by  a  previous  run 
!  You  may  wish  to  delete  this  information  before  starting  a  run 
!  where  you  will  (interactively)  specify  a  different  iteration 
!  mode.  Edit  this  table  only  if  you  really  know  your  model! 
INVARS        504050607  27  5 
TARGS         5  29  6  31    1  32   1  33   1  34   1 
SPECIALS     0 

PLTVAR       702-1    1010203040571  20   1 
!  Plot  start,  end,  and  step  values.  May  be  edited  if  you  wish. 
!  Outer  Loop:  I  Inner  Loop  . 

410.0  427.0         0.2000  1.000  1.000  1.000 
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APPENDIX  A.2  THE  DELTAE  INPUT  FILE  FOR  THE  HIGH  FREQUENCY  MODE 
OF  THE  ANNULAR  THERMO  ACOUSTIC  PRIME  MOVER 


TITLE       Annular  resonator  (file  ::  to2dTfhighC.in) 

!Created@  12:00:08  02-MAY-97  with  DeltaE  Vers.  3.0bl  for  the  Power  Macintosh 


BEGIN       Initial 

1.0130E+05 

a  Mean  P 

Pa 

0.0000 

Alpl@0  G(0d) 

P 

430.5 

b  Freq. 

Hz 

0.0000 

B  Ph(p)0  G(  Oe) 

P 

293.0 

c  T-beg 

K 

0.0000 

C  IUI@0  G(  Of) 

P 

500.0 

d  lpl@0 

Pa      G 

0.0000 

D  Ph(U)0  G(  Og) 

P 

0.0000 

e  Ph(p)0 

deg     G 

0.0000 

E  Heatln  G(27e) 

P 

2.9000E-04 

f  IUI@0 

mA3/s    G 

0.0000 

g  Ph(U)0 

deg     G 

air 

Gas  type 

ideal 

i 

Solid  type 

1   

SOFTEND 

1 

0.0000 

a  Re(Z) 

(t) 

0.0000 

A  Ipl         Pa 

0.0000 

b  Im(Z) 

(0 

0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 

B  Ph(p)       deg 
C  IUI       mA3/s 
D  Ph(U)       deg 
E  Hdot       W 
F  Work       W 
G  Re(Z) 

sameas  0  Gas  type 

0.0000 

H  Im(Z) 

ideal       Solid  type 

—    2  — 

0.0000 

I    T         K 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

A  Ipl         Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p)       deg 

3.1788E-02 

c  Length 

m 

0.0000 
0.0000 

C  IUI       mA3/s 
D  Ph(U)       deg 

sameas  0 

Gas  type 

0.0000 

E  Hdot       W 

ideal 

Solid  type 

0.0000 

F  Work       W 
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! - 3 

ISODUCT     Duct 
2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

! - 4 

ISODUCT     Duct 
2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

! - 5 

ISODUCT     Duct 
2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

! — 6 

ISODUCT     Duct 
2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

! —    7 

ISODUCT     Duct 
2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 


0.0000 

Alpl 

Pa 

0.0000 

B  Ph(p) 

deg 

0.0000 

C  IUI      mA3/s 

0.0000 

D  Ph(U) 

deg 

0.0000 

EHdot 

W 

0.0000 

FWork 

W 

0.0000 

Alpl 

Pa 

0.0000 

B  Ph(p) 

deg 

0.0000 

CIUI      i 

■nA3/s 

0.0000 

D  Ph(U) 

deg 

0.0000 

EHdot 

W 

0.0000 

FWork 

W 

0.0000 

Alpl 

Pa 

0.0000 

B  Ph(p) 

deg 

0.0000 

CIUI      i 

mA3/s 

0.0000 

D  Ph(U) 

deg 

0.0000 

EHdot 

W 

0.0000 

FWork 

W 

0.0000 

Alpl 

Pa 

0.0000 

B  Ph(p) 

deg 

0.0000 

CIUI      i 

mA3/s 

0.0000 

D  Ph(U) 

deg 

0.0000 

EHdot 

W 

0.0000 

FWork 

W 

0.0000 

Alpl 

Pa 

0.0000 

B  Ph(p) 

deg 
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3.1788E-02 

c  Length 

m 

0.0000 

CIUI       i 

rnA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

i 

Solid  type 

...     s 

0.0000 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI       i 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

1 

Solid  type 

...    9  ... 

0.0000 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

! 

Solid  type 

—   10  — 

0.0000 

FWork 

W 

ISODUCT 

Duct 

i  \j 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI      i 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

t 

Solid  type 

1 1   __ 

0.0000 

FWork 

w 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

Solid  type 

0.0000 

FWork 

W 
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ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI       i 

nA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

i 

Solid  type 

1  T,    

0.0000 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI      i 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

i _ 

Solid  type 

—  14  .... 

0.0000 

FWork 

W 

VDUCER 

Driver 

1.0000E-09 

a  Re(Ze) 

ohms 

0.0000 

Alpl 

Pa 

0.0000 

b  Im(Ze) 

0.0000 

B  Ph(p) 

deg 

1 .0000E+04 

c  Re(Tl)  V-s/mA3 

0.0000 

CIUI      i 

mA3/s 

0.0000 

d  Im(Tl)  V-s/mA3 

0.0000 

D  Ph(U) 

deg 

-1.0000E+04 

■     e  Re(T2) 

Pa/A 

0.0000 

EHdot 

W 

0.0000 

f  Im(T2) 

Pa/A 

0.0000 

FWork 

W 

1.0000E-09 

g  Re(Zm)  Pa-s/mA3 

0.0000 

G  Workln     W 

1 .0000E-09 

h  Im(Zm)  Pa-s/mA3 

0.0000 

H  Volts 

V 

1.000 

i  AplVol 

V 

0.0000 

I  Amps 

V 

0.0000 

J  Ph(Ze) 

deg 

sameas  0 

Gas  type 

0.0000 

KlUxl 

mA3/s 

ideal 

Solid  type 

0.0000 

L  Ph(-Ux    deg 

i 

1  s   ... 

ISODUCT 

Constriction 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI 

mA3/s 

0.0000 

D  Ph(U) 

deg 
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sameas  0  Gas  type 

ideal  Solid  type 

! - 1 6 

ISODUCT     Duct  7 

2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

! 1 7 

ISODUCT     Duct  8 

2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

! 1 8 

ISODUCT     Duct 
2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 

sameas  0  Gas  type 

ideal  Solid  type 

! 1 9 

ISODUCT     Duct 
2.6300E-03       a  Area       mA2 
0.2050  b  Perim       m 

3.1788E-02       c  Length      m 


0.0000 
0.0000 


EHdot 
FWork 


W 
W 


0.0000  A  Ipl         Pa 

0.0000  B  Ph(p)       deg 

0.0000  C  IUI       mA3/s 

0.0000  D  Ph(U)       deg 

0.0000  E  Hdot       W 

0.0000  F  Work       W 


0.0000  A  Ipl         Pa 

0.0000  B  Ph(p)       deg 

0.0000  C  IUI       mA3/s 

0.0000  D  Ph(U)       deg 

0.0000  E  Hdot       W 

0.0000  F  Work       W 


Gas  type 
Solid  type 


sameas  0 
ideal 

! _. __ 

ISODUCT     Constriction 


0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 


A  Ipl         Pa 
B  Ph(p)       deg 
C  IUI       mA3/s 
D  Ph(U)       deg 
E  Hdot       W 
F  Work       W 


0.0000  A  Ipl        Pa 

0.0000  B  Ph(p)       deg 

0.0000  C  IUI       mA3/s 

0.0000  D  Ph(U)       deg 

0.0000  E  Hdot       W 

0.0000  F  Work       W 


20 
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2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

i 

Solid  type 

.—  21  - 

0.0000 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

i 

Solid  type 

.—  22  -- 

0.0000 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

? 

Solid  type 

....  23  -- 

0.0000 

FWork 

W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

Alpl 

Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p) 

deg 

3.1788E-02 

c  Length 

m 

0.0000 

CIUI 

mA3/s 

0.0000 

D  Ph(U) 

deg 

sameas  0 

Gas  type 

0.0000 

EHdot 

W 

ideal 

i 

Solid  type 

. 94    .. 

0.0000 

FWork 

W 

ISODUCT     Duct 

2.6300E-03       a  Area  mA2 

0.2050               b  Perim  m 

3.1788E-02       c  Length  m 


0.0000  A  Ipl         Pa 

0.0000  B  Ph(p)       deg 

0.0000  C  IUI       mA3/s 

0.0000  D  Ph(U)       deg 
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sameas  0 

Gas  type 

0.0000 

E  Hdot       W 

ideal 

i 

Solid  type 

—  25  — 

0.0000 

F  Work       W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

A  Ipl        Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p)       deg 

3.1788E-02 

c  Length 

m 

0.0000 
0.0000 

C  IUI       mA3/s 
D  Ph(U)       deg 

sameas  0 

Gas  type 

0.0000 

E  Hdot       W 

ideal 

t_ _ 

Solid  type 

—  26  — 

0.0000 

F  Work       W 

ISODUCT 

Duct 

2.6300E-03 

a  Area 

mA2 

0.0000 

A  Ipl        Pa 

0.2050 

b  Perim 

m 

0.0000 

B  Ph(p)       deg 

3.1788E-02 

c  Length 

m 

0.0000 
0.0000 

C  IUI       mA3/s 
D  Ph(U)       deg 

sameas  0 

Gas  type 

0.0000 

E  Hdot       W 

ideal 

i 

Solid  type 

....  27  — 

0.0000 

F  Work       W 

HXFRST 

Ambient  HX 

sameas  3a 

a  Area 

mA2 

23.02 

A  Ipl         Pa 

0.6290 

b  GasA/A 

103.2 

B  Ph(p) 

deg 

5.0800E-03 

c  Length 

m 

6.6930E-04       C IUI      m 

A3/s 

1.7018E-03 

dyO 

m 

38.94 

D  Ph(U) 

deg 

-1.8309E-03 

e  Heatln 

W      G 

1 .5886E-03       E  Hdot 

W 

293.0 

fEst-T 

K      (t) 

3.3455E-03       FWork 

W 

sameas  0 

Gas  type 

-1.8309E-03      GHeat 

w 

copper 

i 

Solid  type 

.—  28  — 

293.0 

H  MetalT 

K 

STKSLAB 

Prime  Mover  Stack 

2.6300E-03 

a  Area 

mA2 

16.26 

A  Ipl         Pa 

0.6220 

b  GasA/A 

-43.67 

B  Ph(p) 

deg 

2.4000E-02 

c  Length 

m 

6.7080E-04       C  IUI      rr 

iA3/s 

3.0480E-04 

dyO 

m 

38.51 

D  Ph(U) 

deg 

7.5000E-05 

e  Lplate 

m 

1.5886E-03        EHdot 
7.4220E-04       FWork 

w 

W 
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293.0 

G  T-beg        K 

sameas  0 

Gas  type 

293.0 

H  T-end        K 

kapton 

i 

Solid  type 

29  — -- 

-2.6033 

E-03  I  StkWrk 

HXLAST      Cold  HX 

sameas  3a 

a  Area       mA^ 

> 

16.63 

A  Ipl         Pa 

0.7050 

b  GasA/A 

-43.90 

B  Ph(p)       deg 

3.0480E-04 

c  Length      m 

6.7062E-04       C  IUI      mA3/s 

1.2540E-03 

d  yO          m 

38.51 

D  Ph(U)       deg 

0.0000 

e  Heatln      W 

(t) 

7.3668E-04       E  Hdot       W 

293.0 

fEst-T       K 

=29H? 

7.3668E-04       FWork      W 

sameas  0 

Gas  type 

-8.5I95E-04        GHeat        W 

copper 

i 

Solid  type 

30  ----- 

293.0 

H  MetalT      K 

SOFTEND 

0.0000 

a  Re(Z) 

(0 

0.0000 

A  Ipl         Pa 

0.0000 

b  Im(Z) 

(t) 

0.0000 
0.0000 
0.0000 
0.0000 
0.0000 
0.0000 

B  Ph(p)       deg 
C  IUI       mA3/s 
D  Ph(U)       deg 
E  Hdot       W 
F  Work       W 
G  Re(Z) 

sameas  0 

Gas  type 

0.0000 

H  Im(Z) 

ideal 

! 

Solid  type 

^1 

0.0000 

I    T        K 

DIFFTAR 

Ipl    mismatch 

0.0000 

a  TargDi 

=31A? 

8.0065E-05     A  D1-D2 

1A  b  DIAddr 

30A  c  D2Addr 

i 

32  — - 

DIFFTAR 

Ph(p)  mismatch 

0.0000 

a  TargDi 

=32A?  3.3780E-04      A  D1-D2 

IB  b  DIAddr 

30B  c  D2Addr 

i 

33 

D1FF1AR 

IUI    mismatch 

w 
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0.0000  aTargDi  =33 A?      1.3884E-09    A  D1-D2 

1C  b  DIAddr 

30C  c  D2Addr 

! 34  ----- 

DIFFTAR  Ph(U)  mismatch 

0.0000  aTargDi  =34A?     -1.0691E-04  A  D1-D2 

ID  b  DIAddr 

30D  c  D2Addr 

!  The  restart  information  below  was  generated  by  a  previous  run 
!  You  may  wish  to  delete  this  information  before  starting  a  run 
!  where  you  will  (interactively)  specify  a  different  iteration 
!  mode.  Edit  this  table  only  if  you  really  know  your  model! 
INVARS        504050607  27  5 
TARGS         5  29  6  31    1  32   1  33   1  34   1 
SPECIALS     0 
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APPENDIX  B.  THE  MATHEMATICA  PROGRAM  FOR  THE  LEAST 
SQUARE  FIT  TO  THE  FREQUENCY  RESPONSE 


(*  This  mathematica  program  read  the  plotting  data  file  generated  by  *) 
(*  DeltaE  and  perform  a  least  square  fit  to  find  the  resonance  frequency        *) 
Clear  Ali["GlobaT*"] 

(*  Read  the  data  file  generated  by  the  DeltaE  program.  *) 
(*  In  this  example,  the  data  file  is  flowDT100.de  *) 

datfile=ReadList["flowDT100.de",Table[Number,{8 }]];(*  There  are  eight  columns  in  *) 

(*  the  data  file  *) 

(*  Select  the  data  range  used  in  the  least  square  fit,  column  1  is  the  frequencies  *) 
(*  and  column  7  is  the  pressure  amplitude  *) 

data=datfile[[Range[2, 19], { 1 ,7 }]]  (*  Determine  the  bandwidth  of  fitting        *) 
(*  Initial  guesses  of  A,f0,  and  Q  *) 
soln  =  {A->150,f0->437,Q->100}; 

Attributes[y]= { Listable } ; 

y[f_,A_,4_,QJ  :=  A  (f/f0)*2  Abs[Cot[Pi  (f/f0-I/(2  Q))]/(Pi  (f/f0-I/(2  Q)))]; 

Attributes  [yy]  =  {Listable}; 

yy[f J  =  y[f,A,f0,Q]  /.soln; 

xlist  =  Transpose[data][[l]]; 

ylist  =  Transpose[data][[2]]; 

rmserrortA^f^QJ  :=  Sqrt[  Apply [Plus,(ylist  -  y [xlist,  A,f0,Q])A2]]; 

rmserror[A,f0,Q]/.soln; 

FindMinimum[rmserror[A,f0,Q],{A,{  100,200} },{ f0,  {435,439  }},{Q,{  90,1 10}}] 

(*  A  is  an  arbitrary  constant,  f0  the  resonance  frequency,  Q  the  quality  factor  *) 

(*  End  of  program  *) 
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APPENDIX  C.  THE  MATLAB  PROGRAM  FOR  THE  ANNULAR 
THERMOACOUSTIC  PRIME  MOVER 

Function  updateF.m:  The  function  performs  the  actual  iteration  to  find  the  solutions  for 
the  constricted  annular  prime  mover. 

%  Begin  of  function  % 

function    [y,iterations]=updateF(y); 

% 

%  usage  [newy, iterations]  =  updateF(y) 

%  This  function  implements  picards  method  for  updating  y 

%  solving  f(y,x)=0  for  an  annular  resonator 

%  The  input  y  vector  is  a  8  by  1  vector:  [Tm;pre;pim;Ure;Uim;fre;fim;H2dot] 

%  only  the  last  5  elements  (i.e.  Ure,  Uim,  fre,  fim  and  H2)  are  updated 

%  newy  =  y  -  inv(J)*f(y,x); 

%  where  J  is  the  Jacobian  of  f(y) 

%  REQUIRES:  fdjacF.m,  funcvF.m,  computeL.m, 

%  Note:  Before  run  this  program,  make  sure  to  adjust  the  area  ratio  and  select  the 

%  right  temperature  (Just  active  the  line  for  the  temperature  profile) 

%  Declare  the  global  variables 

global  T2  T3  T4  T5  T6  T7  T8  T9  T10  Tl  1  Th  Tc  Tmatrix  RsLa  Las  R  Rst  r 

r  =  2.565e-2;  %  Equivalent  radius  of  the  annular  resonator,  m 

Rs  =  0.7;  %  The  area  ratio  for  the  constriction 

%  Compute  the  equivalent  inductance  of  constriction 

d  =  0.0529;  %  Width  of  annulus,  m 

h  =  0.0497;  %  Height  of  annulus,  m 

dc  =  d*sqrt(Rs);  %  Width  of  the  constriction,  m 

he  =  h*sqrt(Rs);  %  Height  of  annulus,  m 

dst  =  4.28e-2;  %  Width  of  a  slit  in  the  stack  m 

hst  =  6.12e-4;  %  Height  of  a  slit  in  the  stack  m 

La  =  computeL(d,h,dc,hc);        %  Compute  the  equivalent  inductance  of  constriction 

Las  =  computeL(d,h,dst,hst);      %  Compute  the  equivalent  inductance  of  a  slit  in  the 

%  stack 
R  =  sqrt((dc*hc)./(d.*h)); 
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%  Heater  voltage 
%  heater  =  OV 
%  heater  =  8V 
%  heater  =  11V 


%  heater  =  OV 
%  heater  =  8V 
%  heater  =  11V, 
%  heater  =  14V, 


Rst  =  sqrt((dst*hst)./(d.*h));      %  Compute  the  equivalent  area  ratio  for  one  slit 
%  The  measured  temperature  profiles  for  all  the  measurements 
%  When  run  this  main  program,  select  the  profile  based  on  the  porosity 
%  The  low  mode  and  high  mode  share  the  same  temperature  profile 
%  These  temperatures  are  in  °C 
%  Low  mode  and  high  mode  with  constriction  Rs=0.1 
%  Tmatrix=[21.0  21  21  21  21  21  21  21  21  20.6]; 
%  Tmatrix=[66.9  23.4  22.5  58.8  22.2  21.9  22.2  32.5  25.7  21.4]; 
%  Tmatrix=[  100.7  25.8  24.3  86.8  23.8  23.2  23.7  41.2  30.4  22.6]; 
%  Tmatrix=[  137.8  28.5  26.2  119.2  25.4  24.4  25.2  51.7  36.2  23.6];    %  heater  =  14V 
%  Tmatrix=[218.6  36.6  32.4  194.5  30.6  28.7  30.2  82.1  51.8  27.3];    %  heater  =  20V 
%  Tmatrix=[278.6  43.9  38  253.1  35.2  32.3  34.5  113.1  66.2  30.2];     %  heater  =  25V 
%  Low  Mode  and  high  mode  with  constriction  Rs=0.3 
%  Tmatrix=[29.7  29.5  29.4  29.6  29.5  29.4  29.4  29.3  29.5  29.5]; 
%  Tmatrix=[74.3  30  29.1  66.8  28.9  28.5  28.8  37.8  32.2  28.4]; 
%  Tmatrix=[106.6  31.1  29.6  93.8  29.1  28.4  29.0  44.5  35.1  28.2]; 
%  Tmatrix=[  143.9  33.1  30.7  126  29.9  28.9  29.7  54.2  39.4  28.4]; 
%  Tmatrix=[220.7  38.5  34.4  196.2  32.6  30.6  32.1  81.4  50.8  29.5];    %  heater  =  20V 
%  Tmatrix=[278.4  44.9  39.1  251.5  36.3  33.2  35.5  110.3  63.1  31.5];%  heater  =  25V 
%  Low  and  high  Mode  with  constriction  Rs=0.7 

%  Tmatrix=[23.5  23.5    23.5  23.5  23.5  23.5  23.5  23.5  23.4];  %  heater  =  0V 

%  Tmatrix=[74.9  25.3  24.3  65.7  24  23.6  23.9  55.1  23.7  23.4];  %  heater  =  8V 

%  Tmatrix=[108.5  27.2  25.5  94  24.9  24.2  24.8  74.5  30.4  23.9];  %  heater  =  11V 
%  Tmatrix=[147.7  29.8  27.3  128.2  26.3  25.1  26.0  100.1  34.5  24.6];%  heater  =  14V 
%  Tmatrix=[221.9  35.6  31.2  196  29.4  27  28.8  157.1  43.6  26.2];  %  heater  =  20V 
%  Tmatrix=[280.7  42.5  36.2  252.9  33.3  29.7  32.4  219.2  53.5  29.5];%  heater  =  25V 
Tmatrix  =  Tmatrix+273;  %  Change  temperature  to  Kelvin 

T2=Tmatrix(l); 
T3=Tmatrix(2); 
T4=Tmatrix(3); 
T5=Tmatrix(4); 
T6=Tmatrix(5); 
T7=Tmatrix(6); 
T8=Tmatrix(7); 
T9=Tmatrix(8); 
T10=Tmatrix(9); 
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Tll=Tmatrix(10); 
Th  =  (2*T2+T5)/3; 
Tc  =  (2*T3+T4)/3; 
Nmax  =  100; 

y=y(0; 

y(l,D  =  Th; 

iterations=0; 

%  Main  loop  for  iterations 

while(l) 

iterations  =  iterations  +  1 ; 

fvec=funcvF(y); 

cheeky  =  norm(fvec./y(l:5)) 
if  checky<1.0E-3 

fprintf(l,'Done\n'); 

y  =  y' 


%  Weighted  hot  end  targeted  temperature,  K 

%  Weighted  cold  end   temperature,  K 

%  Maximum  number  of  integration 

%  Make  sure  the  input  vector  y  is  a  column  vector 

%  Fixed  the  starting  temperature  at  Th 


%  Compute  the  difference  between  the  targets  and  % 
%  calculated  results  % 

%  Check  the  norm  of  the  ratio  of  the  difference  % 
%  Set  the  converging  threshold  % 


%  Print  out  the  solution  % 

[fvec,yend]  =  funcvF(y); 

yend  =  yend'     %  print  out  the  integration  results  at  the  ends  % 

Q  =  yend(l,6)/(2*yend(l,7))    %  Find  the  quality  factor  % 
elseif  iterations  ==  Nmax  %  Terminate  the  iteration  if  iteration  >  N„ 

fprintf(  1 , 'Exceed  max  #  of  iterations'); 

return 
else 

J  =  fdjacF(y,fvec); 
diffy  =  -J\fvec; 
if  checky<1.0e-l 

diffy  =  diffy/2;  %  Close  to  converge,  adjust  the  step  size  % 
elseif  checky<5.0e-2 

diffy  =  diffy/4; 
elseif  checky<1.0e-2 

diffy  =  diffy/6; 
end 


% 


%  calculate  the  Jacobian  % 
%  Use  the  LU  factorization  % 


newy=y 

newy(4:8)  =  y(4:8)+diffy; 

y=newy; 


%  Perform  a  Newton  iteration  to  update  y  % 
%  only  update  Ure,  Uim,  fre,  fim  and  H2  % 


fprintf(l, 'Computing  Update,  iterations  %4.0f  \n',iterations) 
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end 
end 
% 
%  End  of  updateF.m  % 
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Function  funcvF.m:  The  function  integrate  the  coupled  ordinary  differential  equations  of 
the  annular  prime  mover. 

%  Begin  of  function  % 
function  [fvec,yend]  =  funcvF(yin) 

% 

%  Includes  the  ends  effect  from  high  order  mode  method 

%  This  function  integrates  the  annular  resonator  with  the  stack/heat  exchanger  assembly 

%  at  the  end  of  the  resonator  (including  ambient  and  hot  heat  exchanger) 

%  constriction  locates  between  xconl  to  xcon2 

%  Requires:  ductfund.m,  ductfunc21.m,  ductcon.m, 

%  ,  stackfun.m,  ode45.m,  ambhxfunc.m,  hothxfunc,  zstuff.m 

%  Declare  the  global  variables 

global  T2  T3  T4  T5  T6  T7  T8  T9  T10  Tl  1  Th  Tc  Tmatrix  Rs  La  Las  R  Rst 

yin=yin(:); 

d  =  0.0529;       %  Width  of  annulus  m 

h  =  0.0497;       %  Height  of  annulus  m 

dc  =  d*sqrt(Rs);  %  Width  of  the  constriction  m 

he  =  h*sqrt(Rs);  %  Height  of  constriction  m 

cp=1005;  %  Isobaric  heat  capacity  per  unit  mass 

%  These  are  the  positions  of  the  thermocouples,  m 

xp0=  0; 

xpl  1=0.00293;  %  Thermocouple  #  9 

xp  12=0.084345;  %  Thermocouple  #  10 

xpl 3=0.39615;  %  Thermocouple  #  11 

xpl4=0.707955;  %  Thermocouple  #  7 

xpl5=0.7556;  %  Thermocouple  #  8 

xpl  =0.76292;  %  Ambient  heat  exchanger 

xp2  =  0.768;  %  Prime  mover  stack 

xp3  =  0.7920;  %  Hot  heat  exchanger 

xconl  =0.1 33863;  %  45°  long  constriction,  90°  from  the  center  of  the  stack  assembly 

xcon2=  0.23290; 

xpL=0.7923;  %  Effective  length  of  the  annular  resonator,  m  % 

%  These  are  not  the  resolution  of  the  integration 

%  These  are  the  points  in  which  the  results  are  extracted 
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xnl  =  20;  xn2  =  30; 

xn3  =  50;  xn4  =  10; 

fvec  =  zeros(5,l);  %  Set  up  the  fvec  vector 

xllspan  =  linspace(xpO,xpl  l,xnl); 

xl2span  =  linspace(xpl  I,xpl2,xn2); 

xl3aspan  =  linspace(xpl2,xconl,xn2); 

xl3bspan  =  linspace(xconl,xcon2,xn2); 

xl3cspan  =  linspace(xcon2,xpl3,xn2); 

xl4span  =  Iinspace(xpl3,xpl4,xn3); 

xl5span  =  Iinspace(xpl4,xpl5,xn2); 

xl6span  =  linspace(xpl5,xpl,xnl); 

x2span  =  linspace(xpl,xp2,xnl); 

x3span  =  Iinspace(xp2,xp3,xn2); 

x4span  =  linspace(xp3,xpL,xnl); 

yO  =  yin(l:7);  %  make  sure  yO  is  column  vector  and  only  the  first  7  elements  used  in  a  duct 

yO  =  y0(:);        %  Set  up  the  starting  boundary  conditions 

[xll,yll]  =  ode45('ductfunc2',xllspan,y0);    %  Use  the  MATLAB  built  in  ODE45  function 
yOll  =  yll(xnl,:);        %  Set  up  the  boundary  condition  for  next  integration 
[xl2,yl2]  =  ode45('ductfunc21',xl2span,y011); 

y012  =  yl2(xn2,:);        %  Set  up  the  boundary  condition  for  next  integration 
[xl3a,yl3a]  =  ode45('ductfunc22al',xl3aspan,y012); 
%  From  duct  to  a  constriction,  U  is  continuous  but  p2=p,-U,Za 
gamma  =1.4;  %  Ratio,  isobaric  to  isochoric  specific  heats 

rh  =  1 .2825e-2;  %  Hydraulic  radius  of  the  resonator,  rh=r/2,  m 

rhc  =  rh*sqrt(Rs);  %  Hydraulic  radius  in  the  constriction,  m 

omega  =  2*pi*(yl3a(xn2,6)  +  j.*yl3a(xn2,7));%  Extract  the  complex  frequency 
rho  =  353.065./yl3a(xn2,l);        %  Density  of  air  as  a  function  of  temperature  (K) 
c  =  20.0447. *sqrt(yl3a(xn2,l));  %  Speed  of  sound  in  air  as  a  function  of  temperature  (K) 
mue  =  1.846e-5.*(yl3a(xn2,l)/300).A1.5.*(410.4./(yl3a(xn2,l)+110.4));  %  Viscosity  of  air 
K=2.624e2.*(yl3a(xn2,l)./300).A1.5.*(523.831./(yl3a(xn2,l)+... 

245.4.*exp(27.6./yl3a(xn2,l))));  %  Thermal  conductivity,  Air 

Pr  =  0.60928+0.23017. *exp(-0.0028565.*yl3a(xn2,l));  %  Prandtl  number 
dV  =  sqrt(2.*mue./(rho.*omega));        %  Viscous  penetration  depth  % 
dK  =  sqrt(2.*K./(rho.*cp.*omega));      %  Thermal  penetration  depth  % 
fv  =  (l-j)*dV/(2*rh); 

fvc  =  (l-j)*dV/(2*rhc); 
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fk  =  fv/sqrt(Pr); 

fkc  =  fvc/sqrt(Pr); 

k  =  omega/c  *sqrt((l+(gamma-l)*fk)/(l-fv));      %  Wave  number  in  a  regular  duct 

kc  =  omega/c  *sqrt((l+(gamma-l)*fkc)/(l-fvc));  %  Wave  number  in  the  constriction 

[H,HmO,mymz2,nynz2]=zstuff(d,h,dc,hc); 

Yl=diag(-j*sqrt(mymz2-k.A2)./(k.*rho*c)); 

Y2=diag(-j*sqrt(nynz2-kc.A2)./(kc.*rho*c)); 

Zimag=(HmO.'*((H*Y2*H.'+Yl)\HmO))/(Rs*4*rA2);  %  Compute  the  acoustic  impedance 

%  caused  by  the  constriction,  using 
%  higher  order  mode  method 

Ral  =  rho.*omega.*dV.*(l-R)./R.*(l+(l-R.A2)./(pi*R).*log((l+R)./(l-R)))/(Rs*4*rA2); 

Zal  =  Ral+  Zimag; 

%Zal=0;  %  If  do  not  want  to  included  end  correction,  just  set  the  impedance  to  zero 

%  by  activating  this  line 

Ul  =  yl3a(xn2,4)  +  j.*yl3a(xn2,5);  %  Volume  velocity  at  entrance  of  the  constriction 

y013b  =  yl3a(xn2,:); 

y013b(l,2)  =  y013b(l,2)  -  real(Zal*Ul);         %  Correct  real  part  of  pressure 

y013b(l,3)  =  y013b(l,3)  -  imag(Zal*Ul);      %  Correct  imag  part  of  pressure 

[xl3b,yl3b]  =  ode45('ductcon',xl3bspan,y013b);  %  Constriction  region 

%  From  a  constriction  to  a  duct ,  U  is  continuous  but  p3=p2-U2Za 

rho  =  353.065./yl3b(xn2,l); 

c  =  20.0447. *sqrt(yl3b(xn2,l)); 

mue  =  1.846e-5.*(yl3b(xn2,l)/300).A1.5.*(410.4./(yl3b(xn2,l)+110.4)); 

K  =  2.624e-2.*(yl3b(xn2,l)./300).A1.5.*(523.8306./(yl3b(xn2,l)+... 

245.4. *exp(-.27.6./yl3b(xn2,l))));%  Thermal  conductivity,  Air 

Pr  =  0.60928+0.23017.*exp(-0.0028565.*yl3b(xn2,l));%Prandtl  number 

dV  =  sqrt(2.*mue./(rho.*omega));  %  Viscous  penetration  depth  % 

dK  =  sqrt(2.*K./(rho.*cp.*omega));  %  Thermal  penetration  depth  % 

fv  =  (l-j)*dV/(2*rh); 

fvc  =  (l-j)*dV/(2*rhc); 

fk  =  fv/sqrt(Pr); 

fkc  =  fvc/sqrt(Pr); 

k  =  omega/c*sqrt((l+(gamma-l)*fk)/(l-fv));  %  Wave  number  in  regular  duct 

kc  =  omega/c*sqrt((l+(gamma-l)*fkc)/(l-fvc));  %  Wave  number  in  the  constriction 

Yl=diag(-j*sqrt(mymz2-k.A2)./(k.*rho*c)); 

Y2=diag(-j*sqrt(nynz2-kc.A2)./(kc.*rho*c)); 
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Zimag=(HmO.'*((H*Y2*H.'+Yl)\HmO))/(Rs*4*rA2); 

Ra2  =  rho.*omega.*dV.*(l-R)./R.*(l+(l-R.A2)./(pi*R).*log((l+R)./(l-R)))/(Rs*4*rA2); 

Za2  =  Ra2  +  Zimag; 

%Za2  =  0; 

U2  =  yl3b(xn2,4)  +  j*yl3b(xn2,5);%  Volume  velocity  at  start  of  the  constriction 

y013c  =  yl3b(xn2,:); 

y013c(l,2)  =  y013c(l,2)-real(Za2*U2);  %  Correct  real  part  of  pressure 

y013c(l,3)  =  y013c(l,3)-imag(Za2*U2);         %   Correct  imag  part  of  pressure 

[xl3c,yl3c]  =  Ode45('ductfunc22a2',xl3cspan,y013c); 

y013  =  yl3c(xn2,:); 

[xl4,yl4]  =  Ode45('ductfunc23',xl4span,y013); 

y014  =  yl4(xn3,:); 

[xl5,yl5]  =  ode45('ductfunc24',xl5span,y014); 

y015  =  yl5(xn2,:); 

[xl6,yl6]  =  Ode45('ductfunc25',xl6span,y015); 

y016  =  yl6(xnl,:); 

[x2,y2]  =  ode45('ambhxfunc',x2span,y016);  %  Integrate  in  the  ambient  heat  exchanger 

%  Compute  the  end  corrections  for  the  stack 

%  Activate  these  command  to  include  the  end  corrections  in  the  stack  region 
%nslit  =  57;  %  Number  of  the  slit 

%Aslit  =  2.61936e-5;    %  Area  for  one  slit,  m2 
%y0=3.06e-4;  %  Half  spacing  of  the  stack,  m 

%ds  =  4.28e-2;  %  Width  of  one  slit  in  the  stack  region,  m 

%hs  =  6.12e-4;  %  Height  of  one  slit  in  the  stack  region,  m 

%omega  =  2*pi*(y2(xnl,6)  +  j.*y2(xnl,7));  %  Complex  frequency 
%rho  =  353.065  ./y2(xnl,l); 
%c  =  20.0447.*sqrt(y2(xnl,l)); 

%mue  =  1.846e-5.*(y2(xnl,l)/300).A1.5.*(410.4./(y2(xnl,l)+110.4)); 
%dV  =  sqrt(2.*mue./(rho.*omega));    %  Viscous  penetration  depth  in  regular  duct 
%dK  =  sqrt(2.*K./(rho.*cp.*omega));    %  Thermal  penetration  depth  in  regular  duct 
%fv  =  (l-j)*dV/(2*rh); 
%fk  =  fv/sqrt(Pr); 

%fvs  =  tanh((l+j)*yO/dV)/((l+j)*yO/dV);  %  fv  for  Parallel  plate 

%fks  =  fvs/sqrt(Pr);  %  fk  for  Parallel  plate 

%k  =  omega/c  *sqrt((l+(gamma-l)*fk)/(l-fv));    %  Wave  number  in  regular  duct 
%ks  =  omega/c  *sqrt((l+(gamma-l)*fks)/(l-fvs));%  Wave  number  in  the  constriction 


144 


%[H,HmO,mymz2,nynz2]=zstuff(d,h,ds,hs); 

%Yl=diag(-j*sqrt(mymz2-k.A2)./(k.*rho*c)); 

%Y2=diag(-j*sqrt(nynz2-ks.A2)./(ks.*rho*c)); 

%Zimag=(HmO.'*((H*Y2*H.'+Yl)\HmO))/(Aslit); 

%Rast3  =  rho.*omega.*dV.*(l-Rst)./Rst.*(l+(l-Rst.A2)./(pi*Rst).*log((l+Rst)./(l-... 
Rst)))/(Aslit); 

%Zastl  =  (Rast3  +  Zimag)/nslit;   %  Effective  impedance  for  end  corrections  of  the 

%  stack 

%Ustl  =  y2(xnl,4)  +  j.*y2(xnl,5);  %  Volume  velocity  at  entrance  of  the  constriction 
Zastl=0;  %  No  end  corrections  for  stack  if  this  line  is  active 

y02  =  [y2(xnl,:)  yin(8)];  %  Add  the  8th  element  (i.e.  H2)  in  the  stack  region 

%y02(l,2)  =  y02(l,2)  -  real(Zastl*Ustl);         %  Correct  real  part  of  pressure 

%y02(l,3)  =  y02(l,3)  -  imag(Zastl*Ustl);       %  Correct  imag  part  of  pressure 
[x3,y3]  =  ode45('stackfun',x3span,y02);  %  Integrate  in  the  prime  mover  stack  region 
%  Activate  these  commands  to  include  the  end  corrections  in  the  stack  region 

%Aslit  =  2.61936e-5;  %  Area  for  one  slit 

%rho  =  353.065./y3(xn2,l); 

%c  =  20.0447.*sqrt(y3(xn2,l)); 

%mue  =  1.846e-5.*(y3(xn2,l)/300).A1.5.*(410.4./(y3(xn2,l)+110.4)); 

%dV  =  sqrt(2.*mue./(rho.*omega));    %  Viscous  penetration  depth  in  regular  duct 

%dK  =  sqrt(2.*K./(rho.*cp.*omega));%  Thermal  penetration  depth  in  regular  duct 

%fv  =  (l-j)*dV/(2*rh); 

%fk  =  fv/sqrt(Pr); 

%fvs  =  tanh((l+j)*yO/dV)/((l+j)*yO/dV);         %  fv  for  Parallel  plate 

%fks  =  fvs/sqrt(Pr);  %  fk  for  Parallel  plate 

%k  =  omega/c*sqrt((l+(gamma-l)*fk)/(l-fv));%  Wave  number  in  regular  duct 

%ks  =  omega/c*sqrt((l+(gamma-l)*fks)/(l-fvs));%  Wave  number  in  the  constriction 

%[H,HmO,mymz2,nynz2]=zstuff(d,h,ds,hs); 

%Yl=diag(  -j*sqrt(mymz2-k.A2)./(k.*rho*c)); 

%Y2=diag(  -j*sqrt(nynz2-ks.A2)./(ks.*rho*c)); 

%Zimag=(HmO.'*((H*Y2*H.'+Yl)\HmO))/(Aslit); 

%Rast4  =  rho.*omega.*dV.*(l-Rst)./Rst.*(l+(l-Rst.A2)./(pi*Rst).*log((l+Rst)./(l-... 
Rst)))/(Aslit); 

%Zast2  =  (Rast4  +  Zimag)/nslit;%  Effective  impedance  for  end  corrections  of  stack 
Zast2=0;  %  If  this  line  ic  active,  no  end  corrections  for  the  stack  is  included 

%Ust2  =  y3(xn2,4)  +  j.*y3(xn2,5);       %  Volume  velocity  at  start  of  the  constriction 
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y03  =  y3(xn2,l:7);        rc  only  the  first  7  elements  are  used  in  the  hot  heat  exchanger 
%y03(1.2)  =  y03(1.2i  -  real(Zast2*Ust2);         %   Correct  real  part  of  pressure 
Tcy03(1.3)  =  y03(l,3)  -  imag(Zast2*Ust2);       %   Correct  imag  part  of  pressure 
[\4.\4]  =  ode45('hothxfunc'.x4span,y03); 

vend  =  [y4(,xnl.:)  yin(8)]';        9c  This  is  the  result  at  the  end  of  the  hot  heat  exchanger 
fvec(l.l)  =  yend(l,l)-Th;  *  Find  tne  difference  of  the  targeted  hot  temperature 

fvec(2:5,l)=  yin(2:5,l)-yend(2:5,l);%  Find  the  differences  of  the  target  pre,  pirn,  Ure  ,  Uim 
vend  =  yend(:):  %  Show  results  of  the  integration  at  the  end 

9c 

%   End  of  funcvF.m 
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Function  fdjacF.m:  This  function  computes  the  Jacobian  using  finite  difference. 

%  Begin  of  function  % 

function  df=fdjacF(yin,fin) 

%  USAGE:  df=fdjacF(yin,fin) 

% 

%  Inputs:   y,  funcvF  input  vector  to  compute  Jacobian  about 

%  f  evaluation  of  funcv  at  y 

%  Outputs:  df  =  matrix  with  df(i,j)  =  dfuncv(i),dy(j) 

%  This  function  compute  an  numerical  approximation  to  the  Jacobian  matrix 

%  for  the  function  funcvF 

% 

%  INPUTS:  yin  in  the  original  guess  vector 

%  fin  is  the  original  differences  of  the  target 

%  Requires:  funcvF. m 


yin=yin(:); 

fin=fin(:); 

h=lE-10.*abs(yin)+eps;    %  make  y  perturbation  be  1E-10  of  y's  current  value 

%  Perturb  with  respect  to  Ur,  Us,  fr,  fs  and  H2  .  i.e.  y(4:8) 

%  DON'T  perturb  with  respect  to  Tm,  pr  or  pi 

df=zeros(5,5);  %  Initialize  the  Jacobian  matrix  (5  by  5) 

for  i=4:8 

%  perturb  the  only  ith  element  of  y  % 

ypert=yin; 

ypert(i)=h(i)+yin(i); 

%  find  f  with  the  perturbed  y  % 

pert  =  funcvF(ypert); 

%  Forward  difference  formula  % 

df(:,i-3)  =  (fpert(:)-fin(:))./h(i);  %  Jacobian  array 

fprintf(l,'.'); 

end 

% 

%  End  of  funcvF.m  % 
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Function  ductfunc2.m:  This  function  contains  the  differential  equations  in  a  duct  region. 

%  Begin  of  function  % 

function  ydot  =  ductfunc2(x,y) 

%  This  function  contains  the  differential  equation  in  a  dust  region 

%  dp,/dx  =  -jw*rho*Ul/((l-fv)*Area) 

%  dU,/dx  =  -jw*Area*(l+(gamma-l)*fk)*pl/(rho*cA2),  U,  is  the  volumetric  velocity 

%  The  input  y  is  a  column  vector  with  the  component  of 

%  y=[Tm;pre;pim;Ure;Uim;frc;fim;H2]; 

%  in  a  duct  region,  H2  is  not  used,  only  y(l:7)  is  used  in  duct  region 

% 

% 

global  T2  T3  T4  T5  T6  T7  T8  T9  T10  Tl  1  Th  Tc  r     %  All  input  temperature  in  C 

xpll  =  0.00293;  %  Location  of  the  thermocouple  #9 

yd)  =  Th; 

y=y(l:7); 

r_h  =  r/2  ;  %  Hydraulic  radius  of  the  resonator 

Area  =  4*rA2;  %  Cross-section  area  of  the  duct,  m2 

Tm  =  y(l); 

prl  =  y(2); 

pirn  =  y(3); 

Url  =  y(4); 

Uim  =  y(5); 

frl  =  y(6); 

fim  =  y(7); 

U  =  Url  +  j.*Uim;  %  The  volumetric  velocity,  s/m3 

p  =  prl  +  j.*pim;  %  Acoustic  pressure 

f  =  frl  +  j*fim;  %  Complex  frequency 

w  =  2*pi*f; 

%  Gas  constants  for  air  ,  Tm  in  Kelvin,  at  1  ATM  % 

%  Thermal         conductivity       % 

K  =  2.624e-2*(Tm/300)A1.5*(523.8306/(Tm+245.4*exp(-27.6/Tm))); 

%  Viscosity  of  Air  % 

mue  =   1.846e-5.*(Tm/300).A1.5.*(410.4./(Tm+110.4)); 
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gamma  =1.4;  %  Ratio,  isobaric  to  isochoric  specific  heats 

cp=1005;  %  Isobaric  heat  capacity  per  unit  mass 

%  Compute  gas  properties  and  the  f  function 

[fv,fk,c,rho,K,Ksolid,Prandtl]=airproperty(Tm,w,l); 

Tmdot  =  (T9-Th)/xpll;  %  Set  the  temperature  gradient  along  the  duct 

%  set  up  the  ODE  equation  in  a  duct  region 

tempi  =  -j.*w.*rho.*U./((l-fv)*Area); 

temp2  =  -j.*w.*Area.*(l+(gam-l).*fk).*p./(rho.*c.A2)+... 

((fk-fv)*Tmdot*U)/((  1  -fv)*(  1  -Prandtl)*Tm); 

%  Set  up  the   ordinary  differential  equations 

prdot  =    real(templ); 

pidot  =    imag(templ); 

Urdot  =    real(temp2); 

Uidot  =    imag(temp2); 

frldot  =  0; 

fimdot  =  0; 

ydot  =[Tmdot;prdot;pidot;Urdot;Uidot;frldot;fimdot] ; 

% 

%  End  of  ductfunc2.m  % 
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Function  ambhxfunc.m:  This  function  contains  the  differential  equations  in  the  ambient 
heat  exchanger. 

%  Begin  of  function  % 

function  ydot  =  ambhxfunc(x,y) 

%  This  function  contains  the  differential  equations  in  the  ambient  heat  exchanger 

%  dp,/dx  =  -jw*rho*Ul/((l-fv)*Area) 

%  dU,/dx  =  -jw*Area*(l+(gamma-l)*fk)*pl/(rho*cA2),  U,  is  the  volumetric  velocity 

%  The  input  y  is  a  column  vector  with  the  component  of 

%  y=[Tm;pre;pjm;Ure;U,m;frc;fini;H2]; 

%  in  the  ambient  heat  exchanger  region,  H2  is  not  used 

%  Note  that  only  y(l:7)  is  used  here 


global  T2  T3  T4  T5  T6  T7  T8  T9  T10  Tl  1  Th  Tc 
y  =  y(l:7); 


%  All  input  temperature  in  C 


y=y(0; 

Aambgas  =  1.6547e-03; 

yOamb  =  1.7018e-03; 

y(D=T6; 

Tm  =  y(l); 

prl  =  y(2); 

pirn  =  y(3); 

Url  =  y(4); 

Uim  =  y(5); 

frl  =  y(6); 

fim  =  y(7); 

U  =  Url  +j.*Uim; 

p  =  prl  +  j.*pim; 

f  =  frl  +  j*fim; 

w  =  2*pi*f; 


%  Make  sure  yO  is  a  column  vector 

%  Gas  area,  m2 

%  Half  spacing  of  plates  for  ambient  heat  exchanger,  m 

%  Ambient  heat  HX  at  T6  Kelvin 


%  The  volumetric  velocity,  s/m3 
%  Acoustic  pressure 
%  Complex  frequency 


%  Gas  properties  for  air  ,  Tm  in  Kelvin,  at  1  atm  % 

%  Thermal  conductivity 

K  =  2.624e-2*(Tm/300)A1.5*(523.8306/(Tm+245.4*exp(-27.6/Tm))); 
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c  =  20.0447  *sqrt(Tm);  %  speed  of  sound  in  the  air 

mue  =  1.846e-5.*(Tm/300).A1.5.*(410.4./(Tm+110.4));  %  Viscosity  of  Air 

rho  =  353.065/Tm;        %  Density  of  air 

gam  =1.4;  %  Ratio,  isobaric  to  isochoric  specific  heats 

cp=1005;  %  Isobaric  heat  capacity  per  unit  mass 

deltaK  =  sqrt(2.*K./(rho.*cp.*w));    %  Thermal  penetration  depth,  m 

deltaV  =  sqrt(2.*mue./(rho.*w));  %  Viscous  penetration  depth,  m 

rootPrandtl=sqrt(mue.*cp./K); 

fv  =  tanh((l+j)*yOamb/deltaV)/((l+j)*yOamb/deltaV);%  fv  for  Parallel  plate 

fk  =  fv./rootPrandtl;  %  fk  for  Parallel  plate 

%  set  up  the  ODE  equations  in  the  ambient  heat  exchanger  region 

tempi  =  -j.*w.*rho.*U./((l-fv)*Aambgas); 

temp2  =  -j.*w.*Aambgas.*(l+(gam-l).*fk),*p./(rho.*c.A2); 

Tmdot  =  0;        %  Note  that  it  is  assumed  no  temperature  gradient  across  the  ambient  heat 

%  exchanger 
prdot  =   real(templ); 
pidot  =    imag(templ); 
Urdot  =   real(temp2); 
Uidot  =    imag(temp2); 
frldot  =  0; 
fimdot  =  0; 

ydot  =[Tmdot;prdot;pidot;Urdot;Uidot;frldot;fimdot] ; 
% 
%  End  of  ambhxfunc.m  % 
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Function  ductcon.m:    This  function  contains  the  differential  equations  in  the  constricted 
duct. 

%  Begin  of  function  % 

function  ydot  =  ductcon(x,y) 

%  This  function  contains  the  differential  equations  in  the  constriction 

%  dp,/dx  =  -jw*rho*Ul/((l-fv)*Area) 

%  dU,/dx  =  -jw*Area*(l+(gamma-l)*fk)*pl/(rho*cA2),  U,  is  the  volumetric  velocity 

%  The  input  y  is  a  column  vector  with  the  component  of 

%  y=[Tm;pre;pim;Ure;U,m;fre;fim;H2]; 

%  in  a  duct  region,  H2  is  not  used,  only  y(l:7)  is  used  in  duct  region 


global  T2  T3  T4  T5  T6  T7  T8  T9  T10  Tl  1  Th  Tc  Rs  r 
%  Define  the  locations  of  the  constriction  and  thermoacoustics  (in  m) 
xpl2  =  0.084345;         %  Thermocouple  #  10 

%  Thermocouple  #  1 1 
%  Constriction  locates  at  xconl  to  xcon2 
%  45°  long  constriction,  90°  from  center  of  stack 


xpl3  =  0.39615; 
xconl=0. 133863; 
xcon2=  0.23290; 
y=y(l:7); 
r_h  =  r/2*sqrt(Rs) 
Area  =  Rs*4*rA2; 
Tm  =  y(l); 
prl  =  y(2); 
pirn  =  y(3); 
Url  =  y(4); 
Uim  =  y(5); 
frl  =  y(6); 
fim  =  y(7); 
U  =  Url  +  j.*Uim; 
p  =  prl  +  j.*pim; 
f  =  frl  +  j*fim; 
w  =  2*pi*f; 


%  Hydraulic  radius  of  the  constricted  duct,  m 
%  Cross-section  area  of  the  duct,  m2 


%  Volumetric  velocity,  s/m3 
%  Acoustic  pressure 
%  Complex  frequency 


%  Gas  properties  for  air  ,  Tm  in  Kelvin,  at  1  atm  % 
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cp=1005;  %  Isobaric  heat  capacity  per  unit  mass 

c=20.0447.*sqrt(Tm);  %  Speed  of  sound  in  air 

rho=353.065./Tm;         %  Density  of  air 

gam  =1.4;  %  Ratio,  isobaric  to  isochoric  specific  heats 

%  Thermal  conductivity,  Air,  Tm  in  Kelvin 

K  =  2.624e-2.*(Tm./300).Al  .5.*(523.8306./(Tm+245.4.*exp(-27.6./Tm))); 

mue  =  1.846e-5.*(Tm/300).A1.5.*(410.4./(Tm+110.4));  %  Viscosity  of  Air 

Prandtl  =  0.60928+0.23017. *exp(-0.0028565.*Tm);  %  Prandtl  number 

deltaK  =  sqrt(2.*K./(rho.*cp.*w));  %  Thermal  penetration  depth  in  the  constriction,  m 

deltaV  =  sqrt(2.*mue./(rho.*w));  %  Viscous  penetration  depth  in  the  constriction,  m 

fv  =  (l-j).*deltaV./(2*r_h);  %  fv  of  the  constriction 

fk  =  fv./sqrt(Prandtl);  %  fk  of  the  constriction 

Tmdot  =  (Tll-T10)/(xpl3-xpl2);  %  Temperature  gradient  across  the  constriction 

%  set  up  the  ODE  equation  in  a  constricted  duct  region 

tempi  =  -j.*w.*rho.*U./((l-fv)*Area); 

temp2  =  -j.*w.*Area.*(l+(gam-l).*flc).*p./(rho.*c.A2)+... 

((fk-fv)*Tmdot*U)/((l-fv)*(l-Prandtl)*Tm); 
prdot  =    real(templ); 
pidot  =    imag(templ); 
Urdot  =    real(temp2); 
Uidot  =    imag(temp2); 
frldot  =  0; 
fimdot  =  0; 

ydot  =[Tmdot;prdot;pidot;Urdot;Uidot;frldot;fimdot]; 
%  End  of  ductcon.m  % 
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Function  stackfun.m:     This  function  contains  the  differential  equations  in  the  prime 
mover  stack. 

%  Begin  of  function  % 

function  ydot  =  stackfun(x,y) 

%  This  function  contains  the  differential  equations  in  the  prime  mover  stack 

%  to  integrate  in  the  stack  region 

%  The  input  y  is  a  column  vector  with  the  component  of 

%  y=[Tm;pre;pim;Urc;Uim;frc;fim;H2]; 

% 

%  Requires:  airproperty.m 


global  T2  T3  T4  T5  T6  T7  T8  T9  T10  Til  Th  Tc 

y=y(0; 

Agas  =  1.7995e-3;        %  Total  gas  cross  section  area,  m2,  56  cover  glasses 

AsAgas  =  0.2237;         %  Asolid/Agas  in  the  stack  region 

y(l)  =  Tc;  %  Ambient  side's  temperature  of  the  prime  mover  stack 

Tm  =  y(l); 

prl  =  y(2); 

pim  =  y(3); 

Url  =  y(4); 

Uim  =  y(5); 

frl  =  y(6); 

fim  =  y(7); 

H2  =   y(8);  %  Time-averaged  energy  flux,  constant  along  the  stack 

U  =  Url  +  j.*Uim;  %  Volumetric  velocity,  s/m3 

p  =  prl  +  j.*pim;  %  Acoustic  pressure 

f  =  f rl  +  j*fim;  %  Complex  frequency 

w  =  2*pi*f; 

%  Gas  properties  for  air 

[fv,fk,c,rho,K,Ksolid,Prandtl]  =  airproperty(Tm,w,2);     %  index  ==  2  for  parallel  plates 

%  geometry 
gam  =1.4;  %  Ratio,  isobaric  to  isochoric  specific  heats 

cp=1005;  %  Isobaric  heat  capacity  per  unit  mass 
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%  H2  is  the  time-averaged  energy  flux,  constant  along  the  stack  % 

dTmdxl  =  H2/Agas-0.5/Agas*real(p*U'*(l-((fk-fv,)/((l+Prandtl)*(l-fv'))))); 

dTmdx2  =  0.5*rho*cp*abs(U)A2*imag(fk+Prandtl*fv')/(real(w)*(l-PrandtlA2)*abs(l-.. 

fv)A2*AgasA2)  -  K-  AsAgas*Ksolid; 
dTmdx  =  dTmdxl /dTmdx2;     %  Temperature  gradient  in  the  stack  region 
tempi  =  -j.*w.*rho.*U./((l-fv)*Agas); 
temp2  =  -j.*w.*Agas.*(l+(gam-l).*fk).*p./(rho.*c.A2)+((fk-fv)*dTmdx*U)/((l-fv)*(l- 

Prandtl)*Tm); 
Tmxdot  =  dTmdx; 
prdot  =    real(templ); 
pidot  =    imag(templ); 
Urdot  =   real(temp2); 
Uidot  =    imag(temp2); 
frldot  =  0; 
fimdot  =  0; 

H2dot  =  0;        %  H2  is  constant  through  the  stack 
ydot  =[Tmxdot;prdot;pidot;Urdot;Uidot;frldot;fimdot;H2dot] ; 

%  End  of  stackfun.m  % 
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Function  airproperty.m:  This  function  compute  the  properties  for  air. 

%  Begin  of  function  % 

function   [fv,fk,c,rho,K,Ksolid,Prandtl]  =  airproperty(T,w,index) 

%  This  function  compute  the  air  properties  at  1  atm  and  T  Kelvin 

%  if      index  ==  1  »  BL  approximation  for   regular  duct 

%  index  ==  2  »  Parallel  plate  for  stack  region 

%  index  ==  3  »  BL  approximation  for  constricted  duct 

global  Rs 

r_h  =  2.5654E-2/2;       %  Hydraulic  radius  of  the  resonator  is  r/2 

yO  =  3.06e-4;  %  Half  spacing  of  plates  for  stack,  m 

cp=1005;  %  Isobaric  heat  capacity  per  unit  mass 

c=20.0447.*sqrt(T);     %  Speed  of  sound  in  air 

rho=353.065./T; 

%  Thermal  conductivity,  Air 

K  =  2.624e-2.*(T./300).A1.5.*(523.8306./(T+245.4.*exp(-27.6./T))); 

Ksolid  =  0.2*(l-exp(-T/100));  %  Thermal  conductivity,  kapton 

mue  =  1.846e-5.*(T/300).A1.5.*(410.4./(T+110.4));   %  Viscosity  of  Air 

Prandtl  =  0.60928+0.2301 7. *exp(-0.0028565.*T);     %  Prandtl  number 

deltaK  =  sqrt(2.*K./(rho.*cp.*w));        %  Thermal  penetration  depth,  m 

deltaV  =  sqrt(2.*mue./(rho.*w));  %  Viscous  penetration  depth,  m 

if  index  ==  1 

fv  =  (l-j).*deltaV./(2*r_h);      %  fv  for  boundary  layer  approximation 

elseif  index  ==  2 

fv  =  tanh((l+j)*y0/deltaV)/((l+j)*yO/deltaV);  %  fv  for  Parallel  plate 

elseif  index  ==  3 

fv  =  (l-j).*deltaV./(2*r_h*sqrt(Rs));%  fv  for  boundary  layer  approximation 

%        in  the  constricted  duct 

else 

fprintf(l,'  Wrong  index  number  !!') 
return 
end 
fk  =  fv./sqrt(Prandtl); 

%  End  of  airproperty.m  % 
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Function  plotpUF.m:  This  function  plot  the  acoustic  pressure  distribution  in  the  annular 
prime  mover  prime  mover  stack 

%  Begin  of  function  % 

function  [x,y,pnorm]=  plotpUF(yin) 

%  This  function  plot  the  acoustic  pressure  amplitude  of  the  auunlar  prime  mover 

%  It  uses  the  solution  from  "updataF.m"  and  integrates  "yin"  from  x=0  to  x=xL 

% 

%  Note:     The  first  section  of  this  function  which  is  not  included  here, 

%  contains  the  function  "funcvF.m".    It  uses  funcvF.m  to  integrate 

%  the  solution  from  "updateF.m"  from  xO  to  xL 

% 

%  Put  :funcvF.m  here 

% 

dT  =  Th-Tc;  %  Temperature  difference  across  the  prime  mover  stack,  K 

freq  =  y0(6,l)  +j*y0(7,l);         %  Complex  frequency 

%  Combine  the  solution  vector 

%  x  is  the  position  vector,  y  is  the  solution  matrix 

X=[xll;xl2;xl3a;xl3b;xl3c;xl4;xl5;xl6;x2;x3;x4]; 

y=[yll(:,l:7);yl2(:,l:7);yl3a(:,l:7);yl3b(:,l:7);yl3c(:,l:7);yl4(:,l:7);yl5(:,l:7);... 
yl6(:,l:7);y2(:,l:7);y3(:,l:7);y4(:,l:7)]; 

T  =  y(:,l); 

p  =y(:,2)+j*y(:,3);         %  Complex  acoustic  pressure 

U  =y(:,4)+j*y(:,5);        %  Complex  volume  velocity 

pm=zeros(8,l); 

%  Find  the  theoretical  pressure  amplitude  at  the  locations  of  each  microphones 

%  There  are  8  microphones  used  in  the  measurement 

for  i=l:8 

xm(i)  =  0.03483+(i-l)*9.90375e-2;     %  Microphone  location,  m 
[Y,index]  =  min(abs(x-xm(i)));  %  Compute  the  corresponding  index 

pc(i)=pnorm(index);  %  Compute  the  theoretical  values  at 

%  each  microphones'  locations 

end 
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pc=pc(:); 

%  These  are  the  calibrated  measured  pressure  amplitude  at  each  microphones 

%  For  the  constricted  annular  prime  mover 

%  Note  that  all  these  measurements  are  subjected  to  the  measured  temperature  profile  listed 

%  in  the  function  "updateF.m" 

% 

%  Low  mode,  With  constriction 

%  Rs  =  0.7,  low  mode 

%pm=[34.64;20.63;14.82;31.92;31.18;13.42;12.07;28.41]';  %  Heater  =  0V 

%pm=[35.44;21.3;15;32.63;31.94;13.78;12.16;29.07]';  %  Heater  =  8V 

%pm=[35.72;21.77;14.95;33.05;32.53;14.09;12.17;29.39]';  %  Heater  =  11V 

%pm=[36.26;22.37;15.14;33.77;33.16;14.5;12.3;30.03]';  %  Heater  =14V 

%pm=[37.38;23.76;15.3;35.39;35.24;15.57;12.62;31.45]';  %  Heater  =  20V 

%pm=[37.98;24.63;15.46;36.71;36.71;16.29;13.1;32.65]';  %  Heater  =  25 

%  Rs=0.3,  low  mode 

%pm=[34.89;26.09;12.32;34.67;28.18;10.54;11.6;26.91]';  %  Heater  =  0V 

%pm=[39.74;29.49;14.06;38.56;31.29;11.57;13.08;29.95]';  %  Heater  =V8 

%pm=[43.11;32.27;15.82;42.16;34.18;12.56;14.41;32.77]';  %  Heater  =  11V 

%pm=[47.39;35.68;18.16;46.65;37.88;13.71;16.1;36.31]';  %  Heater  =  14V 

%pm=[58.26;45.57;24.5;60.3;48.71;17.44;21.35;47.11]';  %  Heater  =  20V 

%pm=[55.92;55.06;30.23;74.79;60.53;21.48;27.1;58.8]';  %  Heater  =  25V 

%  Rs=0.1,  low  mode 

%pm=[40.57;37.74;34.13;38.99;28.57;10.43;9.76;26.8]';  %  Heater  =  0V 

%pm=[47.45;44.12;39.87;45.48;33.24;11.91;11.76;31.61]';  %  Heater  =  8V 

%pm=[53.45;50.17;45.87;52.2;38;13.43;13.95;36.4]';  %  Heater  =  11V 

%pm=[62.9;59.89;54.71;62.74;45.34;15.77;17.47;43.7]';  %  Heater  =  14V 

%pm=[91.08;97.23;89.03;104.3;75.09;25.35;32.76;73.7]';  %  Heater  =  20V 

%pm=[117.4;155.5;145.9;178.6;128.7;43.3;60.38;135.1]';  %  Heater  =  25V 

%  High  mode,  With  constriction 
%  Rs=0.7,  high  mode 

%pm=[7.29;12.45;11.79;6.47;4.78;9.75;12.213;6.14]';  %  Heater  =  0V 

pm=[8.03;13.09;11.48;5.6;5.44;10.19;11.67;5.17]';  %  Heater  =  8V 

%pm=[8.14;12.99;11.19;5.43;5.66;10.14;11.3;4.92],;  %  Heater  =  11V 
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%pm=[8.41;13.22;11.19;5.39;5.97;10.3;11.23;4.84]'; 

%pm=[8.7;13.29;10.96;5.38;6.43;10.35;10.94;4.75]'; 

%pm=[8.59;13.07;10.67;5.5;6.65;10.26;10.7;4.79]'; 

%  Rs=0.3,  high  mode 

%pm=[6.86;12.83;6.18;6.49;2.87;9.72;11.99;4.47]'; 

%pm=[7.24;13.41;6.43;6.56;3.05;10;12.12;4.32]'; 

%pm=[7.38;13.64;6.69;6.62;3.19;10.18;12.21;4.22]'; 

%pm=[7.45;13.83;6.95;6.63;3.32;10.34;12.27;4.06]'; 

%pm=[7.41;14.4;7.48;6.81;3.67;10.87;12.71;3.87],; 

%pm=[5.56;14.5;7.61;6.89;3.89;11.16;12.93;3.67]'; 

%  Rs=0.1,  high  mode 

%pm=[8.75;15.87;14.26;8.5;2.63;11.37;12.68;4.99]'; 

%pm=[9;16.27;14.48;8.59;2.83;11.67;12.85;4.76]'; 

%pm=[8.98;16.33;14.65;8.6;2.94;11.79;13;4.56]'; 

%pm=[8.98;16.8;15.07;8.84;3.13;12.21;13.66;4.48]' 

%pm=[7.69;17.31;15.12;8.99;3.5;12.87;15.10;4.17]' 

%pm=[6.39;17.4;14.87;9.46;3.77;13.43;15.79;4.01]' 


%  Heater  =  14V 
%  Heater  =  20V 
%  Heater  =  25V 

%  Heater  =  OV 
%  Heater  =  8V 
%  Heater  =  11V 
%  Heater  =  14V 
%  Heater  =  20V 
%  Heater  =  25V 

%  Heater  =  OV 
%  Heater  =  8V 
%  Heater  =  11V 
%  Heater  =  14V 
%  Heater  =  20V 
%  Heater  =  25V 


%  A  least  fit  is  applied  to  the  predicted  values  and  the  measured  data 

A=pm*pc/(pc'*pc);       %  Find  the  least  square  fit  constant 

pmscal=pm./A;  %  Normalized  microphone  voltage  by  the  constant  A 

%  These  are  the  actual  positions  of  the  microphones,  m 

xpm=[0.03483;0.13387;0.23291;0.33195;0.43099;0.53003;0.62907;0.72811]; 

clg 

lx  =[0.7629  0.7629];  %  Mark  the  starting  location  for  stack  assembly 

ly  =  [0  1.2]; 

lxl=[0.133863  0.133863  ]; 

lyl=[0  1.2]; 

lx2=[0.2329  0.2329  ]; 

ly2=[0  1.2]; 

%  Plot  the  predicted  pressure  amplitude  vs.  the  least  square  of  the  measured  data 

plot(x,abs(p)./max(abs(p)),'b-',xpm,pmscal,'ro',lx,ly,'g~',lxl,lyl,'g:',lx2,ly2,'g:') 

legend('Matlab  program  ','Measured  data','  Stack  region',  'Constriction'); 

%title(['Low  mode  with  end  effect  3  ;  dT  =  ',num2str(dT),'  Kelvin  ;...    %  Title  for  low  mode 

%         freq  =  ',num2str(freq),'  ;  Rs  =  ',num2str(Rs)]) 

title(['High  mode  with  end  effect  3  ;  dT  =  ',num2str(dT),'  Kelvin  ;...      %  Title  for  high  mode 


%  Mark  the  location  of  the  constriction 


%  Mark  the  location  of  the  constriction 
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freq  =  ',num2str(freq),'  ;  Rs  =  ',num2str(Rs)]) 
xlabel('Position,  m') 
ylabel('Normalized  pressure') 
axis([0  0.7923  0  1.2]) 

%  End  of  plotpUF.m  % 
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Function  computeL.m:  This  function  computes  the  equivalent  inductance  for  the  series 
acoustic  impedance  Za(w)  of  a  size  discontinuity  between  a  duct  dxh  and  dcxhc  using 
higher  order  mode  theory. 
This  function  is  contributed  by  Ralph  T.  Muehleisen 

%  Begin  of  function  % 

function  La=computeL(d,h,dc,hc); 

% 

%  Inputs 

%  d  =  depth  of  main  channel 

%  h  =  height  of  main  channel 

%  dc  =  depth  of  constriction 

%  he  =  height  of  constriction 

%  Outputs 

%  La  =  equivalent  inductance  of  constriction 

% 

%  This  function  computes  the  equivalent  inductance  for  the 

%  series  acoustic  impedance  Za(w)  of  a  size  discontinuity  between  a 

%  duct  dxh  and  dcxhc  using  higher  order  mode  theory. 

%  The  discontinuity  is  assumed  to  be  symmetric 

%  in  the  y  direction  and  asymmetric  in  the  z  direction 

%  (i.e.  going  from  alxbl  to  a2xb2  the  open  area  in  duct  2  is 

%   (d-dc)/2<y<(d+dc)/2  and  (h-hc)<z<h 

% 

%  To  use  this  add  the  following  lines  to  your  code 

% 

%  In  the  beginning: 

%  La  =  computeL(d,h,dc,hc) 

%  R=sqrt((dc*hc)./(d.*h)); 

%  Ra  =  rho.*w.*dV.*(l-R)./R.*(  1+  (l-R.A2)./(pi*R).*log((l+R)./(l-R))  ); 

%  Za  =  Ra  +  j*w*La; 

R=sqrt((dc*hc)./(d.*h)); 

%  compute  the  coupling  matrix  H,  using  n2=3  modes  in  small  constriction 
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%  and  select  nl  based  on  relative  sizes 

n2=3; 

nl=round(n2/sqrt(R)); 

%  call  h2d  to  get  the  full  coupling  matrix 
[Hfull,my,mz,ny,nz]  =  h2d(R,nl,n2); 
[MH.NH]  =  size(Hfull); 

H00=  Hfull  (1,1); 

%  extract  out  the  higher  order  mode  terms 
H  =  Hfull(2:MH,2:NH); 

%  Extract  out  the  plane  wave  terms 

HmO  =  Hfull(2:MH,l); 

my  =  my(2:MH); 

mz  =  mz(2:MH); 

ny  =  ny(2:NH); 

nz  =  nz(2:NH); 

mymz2  =  (my.*pi./d).A2  +  (mz.*pi./h).A2; 
nynz2  =  (ny.*pi./dc).A2  +  (nz.*pi./hc).A2; 
%  The  inductance  approximation,  it  assumes  k2  «  my2mz2  or  ny2nz2  and  k~kc 

Yin  =  diag(sqrt(mymz2)); 

Y2n  =  diag(sqrt(nynz2)); 

La  =  rho*  (HmO.1  *inv(H*Y2n*H.'  +  Yln)*HmO)/(hc*dc); 

return 

%  End  of  computeL.m  % 
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Function  zstuff.m:  This  function  computes  the  scattering  matrices  used  in  computation  of 
series  acoustic  impedance  Za(w)  of  a  size  discontinuity  between  a  duct  dxh  and  dcxhc 
using  higher  order  mode  theory. 
This  function  is  contributed  by  Ralph  T.  Muehleisen 

%  Begin  of  function  % 

function  [H,HmO,mymz2,nynz2]=zstuff(d,h,dc,hc); 

%  usage    [H,HmO,mymz2,nynz2]=zstuff(d,h,dc,hc); 

% 

%  Inputs 

%  d  =  depth  of  regular  duct 

%  h  =  height  of  regular  duct 

%  dc  =  depth  of  constriction 

%  he  =  height  of  constriction 

% 

%  Outputs 

%  H  =  higher  order  mode  scattering  matrix 

%  HmO  =  scattering  vector  back  to  plane  wave  modes 

%  mymz2,  nynz2  =  indices  of  higher  order  modes 

% 

%  This  function  computes  the  scattering  matrices  used  in  computation 

%  of  series  acoustic  impedance  Za(w)  of  a  size  discontinuity  between  a 

%  duct  dxh  and  dcxhc  using  higher  order  mode  theory. 

%  The  discontinuity  is  assumed  to  be  symmetric 

%  in  the  y  direction  and  asymmetric  in  the  z  direction 

%  (i.e.  going  from  alxbl  to  a2xb2  the  open  area  in  duct  2  is 

%  d-dc)/2<y<(d+dc)/2  and  (h-hc)<z<h 

% 

%  Add  the  following  lines  to  your  program 

%  somewhere  in  the  beginning  add 

%  [H,HmO,mymz2,nynz2]=zstuff(d,h,dc,hc); 

% 

%  In  the  frequency  dependent  part  add  following.  If  the  constriction  is  on  the  right 

%  Yl=diag(  -j*sqrt(  mymz2-k.A2)./(k*rho*c)); 

%  Y2=diag(  -j*sqrt(nynz2-kc.A2)./(kc.*rho*c)); 
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%  Zimag=(HmO.'*((H*Y2*H.'+Y  1  )\HmO))/Sc; 

%  R=sqrt((dc*hc)./(d.*h)) 

%  Ra  =  rho.*w.*dV.*(l-R)./R.*(l+(l-R.A2)./(pi*R).*log((l+R)./(l-R))); 

%  Za  =  Ra+  Zimag; 


R=sqrt((dc*hc)./(d.*h)); 

%  Compute  the  coupling  matrix  H,  using  n2=3  modes  in  small  constriction 

%  and  select  nl  based  on  relative  sizes 

n2=3; 

nl=round(n2/R); 

if  nl>  10 

nl=10; 
end 

%  Call  h2d.m  to  get  the  full  coupling  matrix 
[Hfull,my,mz,ny,nz]=h2d(R,n  1  ,n2); 
[MH,NH]=size(Hfull); 

H00=Hfull(l,I); 

%  Extract  out  the  higher  order  mode  terms 
H=Hfull(2:MH,2:NH); 

%  extract  out  the  plane  wave  terms 

HmO  =Hfull(2:MH,l); 

my=my(2:MH); 

mz=mz(2:MH); 

ny=ny(2:NH); 

nz=nz(2:NH); 

mymz2  =  (my.*pi./d).A2  +  (mz.*pi./h).A2; 

nynz2=  (ny.*pi./dc).A2  +  (nz.*pi./hc).A2; 

return 

%  End  of  zstuff.m  % 
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Function  h2d.m:  This  function  is  used  to  calculate  the  coupling  matrix  H  for  a 
change  of  size  of  a  square  duct  (axa)->(bxb)  with  R=b/a. 
This  function  is  contributed  by  Ralph  T.  Muehleisen 

%  Begin  of  function  % 

function  [H,my,mz,ny,nz]=h2d(R,Nl,N2) 

% 

%  Inputs: 

%  Ry  =  a2/al  =  width  ratio  of  duct, 

%  Rz  =  b2/bl  =  height  ratio  of  duct 

%  Nl=  #  of  modes  in  large  duct  region  and 

%  N2  =  #  of  modes  in  small  duct  in  one  direction. 

% 

%  Outputs: 

%  H  =  coupling  matrix 

%  Mx  =  large  duct  y  mode  matrix 

%  Mz=  large  duct  z  mode  matrix 

%  Ny  =  small  duct  y  mode  matrix 

%  Ny  =  small  duct  z  mode  matrix 

% 

%  The  routine  uses  the  Nl  and  N2  lowest  frequency  modes  in  each  section. 

%  It  assumes  that  the  y  direction  constriction  is  symmetric  (i.e  b  is 

%  centered  on  a)  and  the  z  constriction  is  asymmetric  (b  is  at  one 

%  edge  of  a 

% 

Nt2=N2.A2; 

Ntl=Nl.A2; 

%  Generate  ny,  nz 

temp=(0:N2-l)'; 

ny=zeros(Nt2,l); 

nz=ny; 

for  i=0:N2-l 

ny(i*N2  +  (temp+l),l)=[i.*ones(size(temp))l; 

nz(i*N2  +  (temp+l),l)=temp; 
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end 


%  now  sort  the  modes  by  frequency 

f2=sqrt(ny.A2  +  nz.A2);  %  Compute  relative  freq  of  each  mode 

[y,I]=sort(f2);  %  Sort  the  freq. 

ny=ny(I);  %  Pull  out  ny  and  nz 

nz=nz(I); 

%  Generate  my,mz 

temp=(0:Nl-l)'; 

my=zeros(Ntl,l); 

mz=my; 

for  i=0:Nl-l 

my(i*Nl  +  (temp+l),l)=[i.*ones(size(temp))]; 

mz(i*Nl  +  (temp+l),l)=temp; 
end 

fl=sqrt(my.A2  +  mz.A2);  %  Compute  relative  freq  of  each  mode 

[y,I]=sort(fl);  %  Sort  the  freq. 

my=my(I);  %  Pull  out  ny  and  nz 

mz=mz(I); 

%  convert  my,mz,ny,nz  into  the  matrices  used  to  compute  H=Hy*Hy; 
[Ny,My]=meshgrid(ny,my); 
[Nz,Mz]=meshgrid(nz,mz); 
Mypio2=My.*pi/2; 
Hy=zeros(size(My)); 
Hz=Hy; 

%  first  compute  H  for  the  symmetric  y  direction 
e=(l-R)/2; 
sme=sin(My*pi*e); 
cme=cos(My*pi*e); 
smepr=sin(My*pi*(e+R)); 
smep2r=sin(My*pi*(e+2*R)); 
%  first  fill  in  the  general  terms 

Hy=-2*R.A(3/2).*My.*(sme-(-l).ANy.*smepr)./(pi*((My.*R).A2-Ny.A2+eps)); 
%  find  the  n=m*R  terms  (including  m=n=0)  and  set  them  to  (-l)A(m+n)*Sqrt(R) 
I=find(abs(Ny-R*My)<10*eps); 
temp=(2*My.*pi.*R.*cme-sme+smep2r)./(4*My*pi*sqrt(R)+eps); 
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Hy(I)=temp(I); 

%  find  the  m=0  terms  and  set  them  to  zero 

I=find(My==0); 

Hy(I)=0*I; 

%  Find  the  n=0  terms,    they  were  computed  above  by  the  general  expression 

%  but  are  sqrt(2)  too  big 

I=find(Ny==0); 

Hy(I)=Hy(I)/sqrt(2); 

%  find  the  m=n=0  terms  and  set  them  equal  to  sqrt(R) 

I=find(My==0  &  Ny==0); 

Hy(I)=sqrt(R)*ones(size(I)); 

%  second  compute  H  for  the  asymmetric  z  direction 

e=(l-R); 

sme=sin(Mz*pi*e); 

cme=cos(Mz*pi*e); 

smepr=sin(Mz*pi*(e+R)); 

smep2r=sin(Mz*pi*(e+2*R)); 

%  first  fill  in  the  general  terms 

Hz=-2*R.A(3/2).*Mz.*(sme-(-l).ANz.*smepr)./(pi*((Mz.*R).A2-Nz.A2+eps)); 

%  find  the  n=m*R  terms  (including  m=n=0)  and  set  them  to  (-l)A(m+n)*Sqrt(R) 

I=find(abs(Nz-R*Mz)<10*eps); 

temp=(2*Mz.*pi.*R.*cme-sme+smep2r)./(4*Mz*pi*sqrt(R)+eps); 

Hz(I)=temp(I); 

%  Find  the  m=0  terms  and  set  them  to  zero 

I=find(Mz==0); 

Hz(I)=0*I; 

%  find  the  n=0  terms,    they  were  computed  above  by  the  general  expression 

%  but  are  sqrt(2)  too  big 

I=find(Nz==0); 

Hz(I)=Hz(I)/sqrt(2); 

%  find  the  m=n=0  terms  and  set  them  equal  to  sqrt(R) 

I=find(Mz==0  &  Nz==0); 

Hz(I)  =  sqrt(R)*ones(size(I)); 

H  =  Hy.*Hz; 

%  End  of  h2d.m  % 
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Function  update2st.m:  This  function  iterates  the  ordinary  differential  equation  for  a  two- 
stack  annular  prime  mover. 

%  Begin  of  function  % 

function    [y,iterations]=update2st(y); 


usage  [newy, iterations]  =  update2st(y) 


%  This  function  implements  picards  method  for  updating  y 

%  solving  f(y,x)=0  for  a  resonator 

%  The  input  y  vector  is  a  8  by  1  vector:  [Tm;pre;p,m;Urc;Uim;fre;f,m;H2a;H2b] 

%  only  the  last  6  elements  (i.e.  Ure,  U,m,  fre,  f,m,  H2a,  and  H2b)  are  updated 

% 

%  newy  =  y  -  inv(J)*f(y,x); 

%  where  J  is  the  Jacobian  of  f(y) 

%  REQUIRES:  fdjac2st.m,  func2st.m 

Nmax  =  50; 

y  =  y(0; 

iterations  =  0; 
while(l) 

iterations  =  iterations  +  1 ; 

fvec=func2st(y); 

cheeky  =  norm(fvec./y(l:6))     %  check  the  norm  of  the  ratio  of  the  diffenence  % 

if  cheeky  <  1.0E-2  %  to  the  solution,  be  sure  to  use  "./"  not  7"  % 

fprintfO/DoneW); 

y  =  y'  %  Print  out  the  solution  % 

[fvec.yend]  =  func2st(y); 

yend  =  yend'  %  print  the  integration  results  at  the  ends  % 

Q  =  yend(l,6)/(2*yend(l,7))    %  Find  the  quality  factor  % 

return 
elseif  iterations  ==  Nmax 

fprintf(  1 , 'Exceed  max  #  of  iterations'); 
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return 
else 

J  =  fdjac2st(y,fvec);       %  calculate  the  Jacobian  % 

diffy  =  -J\fvec;  %  Use  LU  factorization  to  find  the  correction  to  the  guess  % 

if  checky<1.0e-l  %  Reduce  the  step  size  when  close  to  solution  % 

diffy  =  diffy/2; 
elseif  cheeky  <  5.0e-2 
diffy  =  diffy/4; 
elseif  cheeky  <  5.0e-3 
diffy  =  diffy/6; 
end 

newy  =  y; 

newy(4:9)  =  y(4:9)  +  diffy;       %  Do  a  Newton  iteration  to  update  y  % 
y  =  newy  %  only  update  Ure,  Uim,  fre,  fim,  H2a,  and  H2b  % 

fprintf ( 1  ,'Computing  Update,  iterations  %4.0f  \n',iterations) 
end 
end 

%  End  of  update2st.m  % 
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Function  fdjac2st.m:    This  function  computes  the  Jacobian  for  the  two-stack  annular 
prime  mover  using  finite  difference. 

%  Begin  of  function  % 

function  df=fdjac2st(yin,fin) 

%  USAGE:  df=fdjac2st(yin,fin) 

% 

%  INPUTS:   y  funcv  input  vector  to  compute  Jacobian  about 

%  f  evaluation  of  funcv  at  y 

%  OUTPUTS:  df  =  matrix  with  df(i,j)  =  dfuncv(i),dy(j) 

%  This  function  compute  an  numerical  approximation  to  the  Jacobian  matrix 

%  for  the  function  funcv 

%  INPUTS: 

%  yin  in  the  original  guess  vector 

%  fin  is  the  original  differences  of  the  target 

% 

%  Requires:  func2st.m 

yin  =  yin(:); 

fin  =  fin(:); 

h  =  lE-10.*abs(yin)+eps;    %  make  y  perturbation  be  1E-10  of  y's  current  value 

%  Perturb  with  respect  to  Ur,  Ui,  fr,  fi  and  H2a  H2b  .  i.e.  y(4:9) 

%  DON'T  perturb  with  respect  to  Tm,  pr  or  pi 

df  =  zeros(6,6);  %  Initialize  the  Jacobian  matrix  (6  by  6) 

for  I  =  4:9 

%  perturb  the  only  ith  element  of  y 

ypert=yin; 

ypert(i)  =  h(i)+yin(i); 
%  find  f  with  the  perturbed  y 

fpert  =  func2st(ypert); 
%  Forward  difference  formula  % 

df(:,i-3)  =  (fpert(:)-fin(:))./h(i);  %  Jacobian  matrix  % 

fprintf(l,'.'); 
end 
%  End  of  fdjac2st.m  % 
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Function  funcv2st.m:    The  function  integrate  the  coupled  ordinary  differential  equations 
of  the  two-stack  annular  prime  mover. 

%  Begin  of  function  % 

function  [fvec,yend]  =  func2st(yin) 

%  This  function  integrates  coupled  ODE  for  the  two  stack  annular  prime  move 

%  Requires: 

%  duct2stl.m,  hothx2stl.m,  stackl.m,  stack2.m,  ambhx2stl.m,  ambhx2st2.m, 

%  hothx2stl.m,  hothx2st2.m,  ode45.m 

%  Note:  Most  of  these  functions  are  similar  to  the  functions  for  the  constricted 

%  annular  prime  mover.  Therefore,  only  the  main  functions  are  included  here. 


global   Th  Tc  %  All  input  temperature  in  K 

%  Before  to  execute  the  program,  be  sure  to  adjust 

%  the  target  hot  heat  exchanger  temperature  to  the  desired  value 

Th  =  393; 

Tc  =  293; 

yin  =  yin(:); 

xpO  =  0; 

xpl  =  0.16869; 

xp2  =  0.1689948; 

xp3  =  0.1929948; 

xp4  =  0.198075; 

xp5  =  0.76292; 

xp6  =  0.768; 

xp7  =  0.792; 

xpL  =  0.7923; 

xnl  =  50; 

xn2  =  10; 

xn3  =  20; 

fvec  =  zeros(6,l); 

xlspan  =  linspace(xpO,xpl,xnl); 

x2span  =  Iinspace(xpl,xp2,xn2); 

x3span  =  Iinspace(xp2,xp3,xn3); 

x4span  =  Iinspace(xp3,xp4,xn2); 


% 
% 


%  Starting  location  for  the  first  satck  assembly 

%  End  of  the  first  stack  assembly 

%  Starting  location  for  the  second  satck  assembly 

%  End  of  the  second  stack  assembly 
%  Effective  length  of  the  resonator,  m 
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x5span  =  Iinspace(xp4,xp5,xnl); 

x6span  =  Iinspace(xp5,xp6,xn2); 

x7span  =  Iinspace(xp6,xp7,xn3); 

x8span  =  Hnspace(xp7,xpL,xn2); 

% 

%  Make  sure  that  yO  is  column  vector,  only  the  first  7  elements  are  used  in  the  duct% 

y0  =  yin(l:7); 

yO  =  y0(:); 

[xl,yl]  =  ode45('duct2stl',xlspan,y0); 

yOl  =  yl(xnl,:); 

y01(l,l)  =  Th; 

[x2,y2]  =  ode45('hothx2stl',x2span,y01);%  Temperature  at  Th 

y02  =  [y2(xn2,:)  yin(8)]';  %  Add  the  8th  element  (i.e.  H2a)  for  the  ist  stack 

[x3,y3]  =  ode45('stackl',x3span,y02); 

y03  =  y3(xn3,l  :7);  %  only  the  first  7  elements  are  used  in  the  hot  heat  exchanger 

[x4,y4]  =  ode45('ambhx2stl',x4span,y03);%  Target  Tc  at  the  1st  ambient  heat  exchanger 

y04  =  y4(xn2,:); 

[x5,y5]  =  ode45('duct2stl',x5span,y04); 

y05  =  y5(xnl,:); 

y05(l,l)  =  Tc; 

[x6,y6]  =  Ode45('ambhx2st2',x6span,y05);%  Temperature  at  Tc 

y06  =  [y6(xn2,:)  yin(8)  yin(9)]';  %  Add  the  8th  element  (i.e.  H2)  for  the  2nd  stack 

[x7,y7]  =  ode45('stack2',x7span,y06); 

y07  =  y7(xn3,l:7);        %  only  the  first  7  elements  are  used  in  the  hot  heat  exchanger 

[x8,y8]  =  ode45('hothx2st2',x8span,y07);         %  Target  Th  at  2nd  hot  heat  exchanger 

%  Combine  results  of  the  integration 

%  This  is  the  result  at  the  end  of  the  2nd  hot  heat  exchanger 

yend  =  [y8(xn2,:)  yin(8)  yin(9)]'; 

%  Compute  the  target  values  % 

fvec(l,l)  =  y03(l,l)-Tc;  %  Find  the  difference  of  the  targeted  cold  temperature  at  ambhx2stl 

fvec(6,l)  =  yend(l,l)-Th;%  Find  the  difference  of  the  targeted  hot  temperature  at  hothx2st2 

%  Find  the  differences  of  the  target  pre,  pirn,  Ure  and  Uim  % 

fvec(2:5,l)=  yin(2:5,l)-yend(2:5,l); 

yend  =  yend(:);  %  Show  results  of  the  integration  at  the  end 

%  End  of  funcv2st.m  % 
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APPENDIX  D.  PROPERTIES  OF  AF  45  GLASS 

Mechanical  properties: 

Density  at  20  °C  (68  °F):  p  =  2.72  g/cm3 

Young's  modulus:  c  =  66.0  kN/mm2 
Torsional  modulus:  E  =  26.7  kN/mm2 

Poisson's  ratio:  y=  0.235 

Electrical  properties: 

Dielectric  constant  (1  MHz):  £r  =  6.2 

Dielectric  loss  factor  (1  MHz):  tan  8  =  9  x  10"4 

Temperature  for  a  specific  electrical  resistivity  of  108  Q.  cm:  TKl00  =  610  °C  (1 130  °F) 

Thermal  properties: 

Viscosity  and  corresponding  temperatures 


Viscosity 
log  T|  (d  Pas) 

Temperature 
in  °C           (°F) 

Designation 

14.5 

627             (1161) 

Strain  point 

13.0 

663              (1225) 

Annealing  point 

7.60 

883              (1621) 

Softening  point 

Transformation  temperature:  Tg  =  662  °C  (1224  °F) 

Coefficient  of  mean  linear  thermal  expansion:  a20.300  =  4.5  x  10"6  K ' 

Optical  properties: 

Refractive  indices  at  20°C  (68°F):  nc(X  =  546  nm)  =  1.5275,  nd  (X  =  546  nm)  =  1.5255 

Abbe  value:  vo  =  62.2% 

Luminous  transmittance  (glass  thickness  1.1  mm):  xVD65  =  91.2% 
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APPENDIX  E.l  CONSTRUCTION  DRAWINGS  FOR  THE  PRIME  MOVER 

STACK  HOLDER 


Horizontal  bar 


0.065+/-0.002" 


0.300  +/■ 
0.002" 


0.025+/- 
0.002" 
0.065+/- 
0.002" 


2.080+/- 
0.002" 


#51  drill-threaded 

(for  2-56  flat  head  screw) 


0.025"  +/-0.002 


61  slits,  0.020"  deep,  0.008 
wide  with  seperation  of  0.030" 
(center  to  center)  aling  from 
the  very  center  one 


0.040 " 

0.030 " 
0.030 " 


Vertical  bar 


0.113 +/-0.002" 


0.048+/-0.002" 


J 


h—    1.910+/-0.002"- 


0.248  +/- 
0.002" 


thickness:  0.1 00"+0.002" 


7 


0.064+/-0.002" 
0.120+/-0.002" 
0.064+/-0.002" 


#43  drill-through 

(for  2-56  flat  head  screw) 
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APPENDIX  E.2  CONSTRUCTION  DRAWINGS  FOR  THE  AMBIENT 

HEAT  EXCHANGER 


Vertical  bar 


0.048+/-0.002" 


0.200  +/- 
0.002"     — 


1 .960+/ 
0.002" 


0.025+/-0.002" 
0.075+/-0.002" 


drill-threaded 

(for  0-80  flat  head  screw) 


10  slits,  0.025"  deep, 
0.028"  wide  with 
seperation  of  0.160" 
(center  to  center 


0.025"  +/-0.002 


0.160  "(edge to  center) 
0.160  "(center  to  center) 
0.160  "(center  to  center) 


0.095  +/-0.002" 
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APPENDIX  E.2(CONTINUED)  CONSTRUCTION  DRAWINGS  FOR  THE 

AMBIENT  HEAT  EXCHANGER 

Horizontal  bar 


j  j—  2.030+/-0.002" 

0.148  +/-~4£"  T^T 

0.002" 


T 


r 


1.015+/-0.002" 


thickness:  0.1 00"+0.002" 


0.047+/-0.002" 


0.074+/- 
0.002" 


drill-through 
(for  0-80  flat  head  screw) 


Central  bar 


thickness  :  0.200"+0.002" 
0.050+/-0.002" 


1 .758+/-0.002"  j 


0.100"+/-0.002 


drill-threaded  (for  0-80  flat  head  screw) 
both  ends 


10  slits,  0.025"  deep,  0.028"  wide 
with  seperation  of  0.160"  (center  to  center) 
from  edge  to  the  first  slit  is  0.159" 
both  sides 
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APPENDIX  E.3  CONSTRUCTION  DRAWINGS  FOR  THE  HOT  HEAT 


C 


(2)    p 

1.960+/- 
0.002" 


• 


EXCHANGER 

2.080+/-0.002"  j 

1     0.100"+/- 
I     0.002" 


0.130"+/-0.002" 


(1) 


I 
0.130"+/-0.002" 


0.100"+/- 
0.002" 


0.320+/-0.002" 


0.100"+/- 
0.002" 


I 

! 

—\ 
1 

•  •        3 

.  I 

I 

I 
j 

■   :.  I 

;■;;'        ] 

(1).  9  threaded  holes  of  spacing  0.205"(center  to  center)  for  0-80  3/16"  fine  screw 
(2).  8  threaded  holes  of  spacing  0.205"  (center  to  center)  for  0-80  3/16"  fine  screw 


0.263+/-0.002" 


0.205+/-0.002 

T 


0.160+/-0.002" 


0.205+/-0.002 


0.160+/-0.002" 


0.5O0+/-0.0I  0.500+/-0.002 " 

drill  holes  for  00-90  threaded  rod 
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APPENDIX  F.  GRAPHS  OF  MATLAB  AND  MEASURED  RESULTS 

This  section  which  is  divided  into  two  subsections  contains  graphs  of  the  results  of 
the  MATLAB  program  and  measurements.  First,  results  of  the  constricted  annular  prime 
mover  are  presented.  Next,  some  predictied  results  for  the  two  stack  annular  prime  mover 
are  provided. 
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APPENDIX  F.l  RESULTS  OF  THE  CONSTRICTED  ANNULAR  PRIME 

MOVER 


Low  mode,  Resonance  frequency  vs.  DeltaT  ,  Rs  =  0.7 
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-  -  Calculated,  end  corrections  at  constriction  &  stack  (B) 
o     Measured  data 
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Figure  F.l-1  Comparisons  of  the  measured  and  calculated  resonance  frequencies 
of  the  low  mode  for  the  constricted  prime  mover  with  the  porosity  of  0.7. 
Method  A  is  the  conformation  transformation  and  Method  B  is  the  higher  order 
mode. 
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High  mode,  Resonance  frequency  vs.  DeltaT  ,  Rs  =  0.7 
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Figure  F.l-2  Comparisons  of  the  measured  and  calculated  resonance  frequencies 
of  the  high  mode  for  the  constricted  prime  mover  with  the  porosity  of  0.7. 
Method  A  is  the  conformation  transformation  and  Method  B  is  the  higher  order 
mode. 
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Low  mode,  Resonance  frequency  vs.  DeltaT  ,  Rs  =  0.3 
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o     Measured  data 
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Figure  F.l-3  Comparisons  of  the  measured  and  calculated  resonance  frequencies 
of  the  low  mode  for  the  constricted  prime  mover  with  the  porosity  of  0.3. 
Method  A  is  the  conformation  transformation  and  Method  B  is  the  higher  order 
mode. 
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High  mode,  Resonance  frequency  vs.  DeltaT ,  Rs  =  0.3 
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Figure  F.l-4  Comparisons  of  the  measured  and  calculated  resonance  frequencies 
of  the  high  mode  for  the  constricted  prime  mover  with  the  porosity  of  0.3. 
Method  A  is  the  conformation  transformation  and  Method  B  is  the  higher  order 
mode. 
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Low  mode  with  end  effect  2  ;  dT  =  0  Kelvin  ;  freq  =  41 7.9+2. 149i ;  Rs  =  0.7 
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Figure  F.l-5     Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.      The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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Low  mode  with  end  effect  2  ;  dT  =  77.03  Kelvin  ;  f req  =  420.5+1 .801  i ;  Rs  =  0.7 


0.3  0.4  0.5 

Position,  m 


Figure  F.l-6    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=77  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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Low  mode  with  end  effect  2  ;  dT  =  231  Kelvin  ;  freq  =  429.1+0.5389i ;  Rs  =  0.7 
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Figure  F.l-7    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=23 1  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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High  mode  with  end  effect  2  ;  dT  =  0  Kelvin  ;  freq  =  438.7+5.992i ;  Rs  =  0.7 


Matlab  program 

o     Measured  data 
—    Stack  region 
Constriction 


0.3  0.4 

Position,  m 

Figure  F.l-8    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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High  mode  with  end  effect  2  ;  dT  =  77.03  Kelvin  ;  freq  =  445+6.1 53i ;  Rs  =  0.7 

I  I 1 1 1 r 


Matlab  program 

o     Measured  data 
—    Stack  region 
Constriction 


0.3  0.4  0.5 

Position,  m 

Figure  F.l-9    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=77  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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High  mode  with  end  effect  2  ;  dT  =  231  Kelvin  ;  freq  =  460+6.746i ;  Rs  =  0.7 
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Figure  F.l-10    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=231  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.    The  frequency  is 
the  calculated  value. 
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Low  mode  with  end  corrections  ;  dT  =  228.2  Kelvin  ;  freq  =  307.4-1 .165i ;  Rs  =  0.1 

T 


Matlab  program 

o     Measured  data 

Stack  region 

Constriction 
*     Onset  data,  dT=260K 


0.3  0.4  0.5 

Position,  m 

Figure  F.  1-1 1     Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.1)  when  the  driver  is  located  45°  from  the  stack  and  AT=228  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.  The  frequency  is 
the  calculated  value.  Symbol  *  represents  the  mode  shape  above  onset,  the 
measured  frequency  is  312  Hz. 
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Low  mode  with  end  effect  2  ;  dT  =  0.2  Kelvin  ;  freq  =  361 .7+2.61 9i ;  Rs  =  0.3 
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Figure  F.  1-12    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.      The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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Low  mode  with  end  effect  2  ;  dT  =  71 .73  Kelvin  ;  freq  =  362.6+1 .847i ;  Rs  =  0.3 
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Figure  F.l-13    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=72  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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Low  mode  with  end  effect  2  ;  dT  =  226.5  Kelvin  ;  freq  =  369.9+0.2682i ;  Rs  =  0.3 
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Figure  F.l-14    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=227  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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High  mode  with  end  effect  2  ;  dT  =  0.2  Kelvin  ;  freq  =  466.9+6.538i ;  Rs  =  0.3 
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Figure  F.l-15    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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High  mode  with  end  effect  2  ;  dT  =  71 .73  Kelvin  ;  freq  =  468.7+6.1 94i ;  Rs  =  0.3 
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Figure  F.l-16    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=72  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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High  mode  with  end  effect  2  ;  dT  =  226.5  Kelvin  ;  freq  =  479.6+5.239i ;  Rs  =  0.3 
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Figure  F.l-17    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=227  K.    The 

calculated  results  are  based  on  the  higher  order  mode  method.   The  frequency  is 
the  calculated  value. 
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Low  mode  with  end  effect  3  ;  dT  =  0  Kelvin  ;  freq  =  295.6+3.873i ;  Rs  =  0.1 
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Figure  F.l-18    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.1)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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Low  mode  with  end  effect  3  ;  dT  =  228.2  Kelvin  ;  freq  =  306.8-0.6823i ;  Rs  =  0.1 
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Figure  F.  1-19    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.1)  when  the  driver  is  located  45°  from  the  stack  and  AT=228  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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High  mode  with  end  effect  3  ; 

dT  =  0  Kelvin  ;  freq  =  468.8+9.638i ;  Rs  =  0.1 
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Figure  F.l-20    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.1)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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High  mode  with  end  effect  3  ;  dT  =  228.2  Kelvin  ;  freq  =  488.9+8.246i ;  Rs  =  0.1 
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Figure  F.l-21     Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.1)  when  the  driver  is  located  45°  from  the  stack  and  AT=228  K.    The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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Low  mode  with  end  effect  3  ;  dT  =  0.2  Kelvin  ;  freq  =  361 .3+2.881  i ;  Rs  =  0.3 
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Figure  F.l-22    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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Low  mode  with  end  effect  3  ;  dT  =  226.5  Kelvin  ;  freq  =  369.6+0.4252i ;  Rs  =  0.3 
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Figure  F.  1  -23    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=227  K.    The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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High  mode  with  end  effect  3  ;  dT  =  0.2  Kelvin  ;  freq  =  461.9+9.625i ;  Rs  =  0.3 
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Figure  F.l-24    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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High  mode  with  end  effect  3  ;  dT  =  226.5  Kelvin  ;  freq  =  475.5+8.55i ;  Rs  =  0.3 
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Figure  F.  1-25    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.3)  when  the  driver  is  located  45°  from  the  stack  and  AT=227  K.    The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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Low  mode  with  end  effect  3  ;  dT  =  0  Kelvin  ;  freq  =  41 7.8+2.21 7i ;  Rs  =  0.7 
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Figure  F.l-26    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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Low  mode  with  end  effect  3  ;  dT  =  231  Kelvin  ;  freq  =  428.8+0.4826i ;  Rs  =  0.7 
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Figure  F.l-27    Mode  shape  of  the  low  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=231  K.    The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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High  mode  with  end  effect  3  ;  dT  =  0  Kelvin  ;  freq  =  434.2+8.779i ;  Rs  =  0.7 
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Figure  F.l-28    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=0  K.     The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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High  mode  with  end  effect  3  ;  dT  =  231  Kelvin  ;  freq  =  456.3+9.995i ;  Rs  =  0.7 
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Figure  F.l-29    Mode  shape  of  the  high  mode  of  the  constricted  prime  mover 

(Rs=0.7)  when  the  driver  is  located  45°  from  the  stack  and  AT=231  K.    The 

calculated  results  are  based  on  the  higher  order  mode  method  and  the  end 
corrections  of  the  stack  are  included.  The  frequency  is  the  calculated  value. 
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APPENDIX  F.2  RESULTS  OF  THE  TWO  STACK  ANNULAR 

PRIME  MOVER 


Mode  shape  for  of  the  two  stack  prime  mover,  dT  =  25K;  freq  =  425+7.421  i 
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Figure  F.2-1  The  calculated  mode  shape  of  the  high  mode  of  the  two  stack 
annular  prime  mover  at  AT  =  25  K.  Case  1:  The  duct  between  the  two  hot  heat 
exchangers  is  held  at  Ta. 
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Mode  shape  for  of  the  two  stack  prime  mover,  dT  =  -5.466e-05K;  freq  =  424+7.234i 
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Figure  F.2-2  The  calculated  mode  shape  of  the  low  mode  of  the  two  stack 
annular  prime  mover  at  AT  =  0  K.  Case  1:  The  duct  between  the  two  hot  heat 
exchangers  is  held  at  Ta. 
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